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Abstract 

This review covers the main aspects of black hole accretion disk theory. We begin with 
the view that one of the main goals of the theory is to better understand the nature of 
black holes themselves. In this light we discuss how accretion disks might reveal some of the 
unique signatures of strong gravity: the event horizon, the innermost stable circular orbit, 
and the ergosphere. We then review, from a first-principles perspective, the physical processes 
at play in accretion disks. This leads us to the four primary accretion disk models that 
we review: Polish doughnuts (thick disks), Shakura-Sunyaev (thin) disks, slim disks, and 
advection-dominated accretion flows (ADAFs). After presenting the models we discuss issues 
of stability, oscillations, and jets. Following our review of the analytic work, we take a parallel 
approach in reviewing numerical studies of black hole accretion disks. We finish with a few 
select applications that highlight particular astrophysical applications: measurements of black 
hole mass and spin, black hole vs. neutron star accretion disks, black hole accretion disk 
spectral states, and quasi-periodic oscillations (QPOs). 



1 Introduction 



Because of its firm connection to black holes themselves, black hole accretion disk theory belongs to 
the realm of fundamental physics. However, the theory itself employs a complicated maze of fluid- 
dynamics results and several phenomenological estimates and guesses known only to specialists. 
This Living Review aims to give readers a useful guide, a "road map" if you will, through the 
unavoidable complexity of the subject. 

Below we list fourteen "Destinations" on this road map and explain their logical connections. In 
our opinion, these Destinations are the key issues, or landmarks, of the theory being reviewed. The 
particular road map that we present is, of course, biased by our own ideas and research histories. 
However, we are confident that the grand landscape of this field will, nevertheless, shine through 
in the end. 

We start by pointing out that black holes are one of the most remarkable inventions of the 
human mind. Their bizarre properties capture nearly everyone's imagination, from Princeton 
string theorists to Hollywood movie makers. Initially, black holes only existed virtually, as weird 
theoretical and mathematical ideas. Models of them were studied with interest, but their real 
existence was questioned by many, including Albert Einstein himself. This situation changed in 
the latter part of the twentieth century, after unambiguous and robust detections were made of 
several astrophysical black hole "candidates" within our Galaxy and in many others. Thus far, 
twcQ classes of black holes have been observed by astronomers. In our theory-minded Living 
Review we do not give detailed descriptions of their observational properties. Instead, we stress 
their importance by starting our road map from the two classes of observed black holes: 

Destination 1: Quasars and other similar supermassive objects, which are collectively called 
"active galactic nuclei" (or AGN), having masse^] in the range 10 6 M Q < M < 1O 9 M . 
They reside at centers of our and other galaxies. The "active" ones among them are the 
most powerful steady energy sources known in the universe. Many have radiant powers L 
in excess of their corresponding Eddington luminosities)^] The high efficiency of quasars, 
r/ = L/Mc 2 > 0.1, where M is the mass supply rate from accretion ("the accretion rate"), 
is puzzling. Black hole accretion disk theory predicts that L > Lgdd would imply small 
accretion efficiency r\ <C 0.1. However, the famous "Soltan argument," based on quasar 
counts, shows that on a long time average, t ~ ^Hubble, quasars can have both L > L^dd and 
77 ~ 0.1 |290U254j . AGN are described in much greater detail in the authoritative monograph 
by Krolik [163]. 

Destination 2: Microquasars and similar "stellar-mass" black holes, having M ~ 10 M®. The 
term "microquasars" was invented by Mirabel |206j to convey that these objects, in many 
regards, behave like scaled-down versions of quasars. A few tens of them have been found in 
our Galaxy as members of X-ray binaries [184] . The natural scaling (time) ~ (mass) adds 
importance to observations of microquasar variability, because the same processes that takes 
hundreds of years (say) in quasars, takes only minutes in microquasars. Particularly interest- 
ing are spectral state changes, which occur on timescales ~ 1 day (see [258] for a review), and 
quasi-periodic oscillations (QPOs), which have timescales as short as ~ 1 ms 308, 184, 258 . 

Quasars and microquasars are among the most intriguing astrophysical objects ever discovered, and 
the goal of black hole accretion disk theory is, obviously, to explain their observed properties - but 

An interesting, unanswered question is if there are any missing populations of black holes, e.g., "intermediate" 
mass black holes in the range 100 Mq < M < 10 6 M Q [1801 [2031 l90l Wf\ or "primordial" black holes |56l . 

2 Mq = 1.99 X 10 33 [g] denotes the mass of the Sun, used in astrophysics as a mass unit. 

3 Radiant power in astrophysics is traditionally called "luminosity." At the Eddington luminosity, Z/Edd = 
1.2 X 1O 38 A//M0 [erg/s], radiation force balances the gravity of the central object (with mass M). In the case of 
stars, "super-Eddington" luminosities, L > Z/Edd > are n °t possible, as this would mean radiation pressure would 
blow the star apart. 



2 



one hopes for much more! One hopes that observations of quasars and microquasars, together with 
their proper theoretical interpretation, would eventually test the very heart of black hole physics 
itself. When this happens, we may meaningfully constrain our knowledge of the fundamental 
properties of space and timej^] 

To accomplish this goal, black hole accretion disk theory must find ways to filter out from 
the observational data the parts that bear clear black hole signatures from the remaining dirty 
astrophysics parts. The logical structure of our Living Review attempts to reflect as closely as 
possible this main task of the theory. Thus, our review starts (in Section [2]) with the pure gravity 
part, concentrating on the three particular signatures of black hole gravity that are crucial to black 
hole accretion disk theory: 

Destination 3: Event horizon. This is a sphere of radius ~ GM/c 2 surrounding the black hole 
singularity, from within which nothing may emerge — a one-way membrane. Note that this 
means that black holes have no rigid surfaces. This is a unique signature of black holes; other 
relativistic features may be observable around non-black hole objects, specifically sufficiently 
compact neutron stars, but the event horizon is a defining property of black holes. 

Destination 4: Ergosphere. This is a region around a rotating black hole where spacetime itself 
is dragged along in the direction of rotation at a speed greater than the local speed of light 
in relation to the rest of the universe. In this region, negative energy states are possible, 
which means that the rotational energy of the black hole can be tapped through various 
manifestations of the "Penrose process" |242j . 

Destination 5 : Marginally stable orbit (also called the "innermost stable circular orbit" or ISCO). 
This is the smallest circle (r = r ms ) along which free particles may stably orbit around a 
black hole. No stable circular motion is possible for r < r ms . This is a unique feature of 
relativity, as in Newtonian theory, orbits at all radii are possible^] 

As fascinating as these destinations are, they are not easily probed in nature. The radiation we 
receive from quasars and microquasars comes not from the black holes themselves, but instead orig- 
inates in the accretion disks which surround thenj^](see Figure[l]). In these accretion disks, angular 
momentum is gradually removed by some presumably (although not necessarily |48j ) dissipative 
process, causing matter to spiral down into the black hole, converting its gravitational energy into 
heat, and then, by various processes, radiating this energy^] The radiation subsequently leaks 
through the disk, escapes from its surface, and travels along trajectories curved (in space) by the 
strong gravity of the black hole, eventually reaching our telescopes. Bingo! That is what we are 
interested in this Living Review. 

The preceding description of how black hole accretion works is well established in general, yet 
several rather crucial details remain either not sufficiently understood or are too complex to be 
studied even by the most powerful present-day computers. To deal with this seemingly hopeless 
situation, several purely phcnomenological approaches have been adopted. Two of them should 



4 Observations of black holes may eventually cast light on the quantum gravity structure of the physical vac- 
uum [3], possibly constraining string theory |64| and the hypothesis of extra dimensions 21 251, 140], 

5 Strictly speaking, this statement is only true for nearly spherical gravity sources. Higher order (octopolc) 
moments allow for the formation of an ISCO, even in Newtonian theory [15) ■ 

6 Astrophysical black holes do not themselves radiate. The temperature associated with Hawking radiation is 
T H = (hc 3 )/(8nGMk B ). For a stellar -mass black hole Tjj ~ 10 8 [° K]. Thus, Hawking radiation is completely 
suppressed by the thermal bath of the 3 [° K] cosmic background radiation. For supermassive black holes, the 
Hawking temperature is at least five orders of magnitude smaller still. 

7 The rate M* at which matter accretes is therefore regulated by internal torques and radiative processes. Only 
if M t = Mo = const everywhere, with Mo being the outside mass supply, can the accretion process be stationary. 
Since the internally determined M* may change due to instabilities, limit cycles, etc., an occurrence of a really 
long-term steady accretion flow should be considered a fine-tuned eigenstate. Note, too, that in many astrophysical 
situations Mq is also genuinely variable. 
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Figure 1: An artist's rendition of a generic black hole accretion disk and jet. Inset figures include 
a time sequence of radio images from the jet in microquasar, GRS 1915+105 [204] and an optical 
image of the jet in quasar, M87 (Credit: J. A. Biretta et al., Hubble Heritage Team (STScI /AURA), 
NASA). 
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certainly be marked as Destinations, in and of themselves, on the road map, because they are 
among the very basic ingredients of the accretion process: 

Destination 6: Angular momentum transport and energy dissipation. The quasi-steady accretion 
of a particle of mass m through a Keplerian disk from a large outer radius, r out , to an inner 
radius, r; n , requires that the particle give up an amount of energjj^j^ 0.1 mc 2 . To do this, 
the particle must also give up an amount of angular momentum ~ {GMr^t) 1 ^ 2 ■ We will 



see in Section 3.2 that viscous stresses within the fluid can facilitate this (mass transfer in, 
angular momentum transfer out, and energy dissipation). However, the stresses can not come 
from ordinary molecular viscosity, as this is much too weak in astrophysical accretion disks. 
Instead, the stresses likely come from turbulence that acts like an effective viscosity. 

Destination 7: Radiative processes and radiative transfer. These depend on the thermodynamic 
state of matter (electron density, ion density, temperature), its motion, and the magnetic 
field, but most importantly on whether the matter is opaque or transparent to radiation, i.e., 
whether its optical depth is large, t ^> 1, or small, r< 1. In the general case when the matter 
is optically thick r > 1, the accretion disk can be quite luminous and also efficiently cooled 
by radiation. Accretion disks with r < 1 are inefficiently cooled and thus less luminous. 

With these complexities, it is useful to have both analytic theories and numerical simulations 
to facilitate progress in the field. Before getting into details, though, let us simply study the 
parameter space available before us. Let us first consider three generic types of physical processes 
that must occur in black hole accretion disks: "dynamical" processes with a characteristic time- 
scale, tdyn ~ l/fi, where fl is the orbital angular velocity; "thermal" processes with a characteristic 
time-scale, t t h ~ c 2 /^^ 2 , where c s is the sound speed and is the kinematic viscosity; and 
"viscous" processes with a characteristic time-scale, i v is ~ r 2 where r is the radial distance 
from the black hole. One excepts that, typically^] 

tdyn "C Uh "C ivis- (1) 

Because dynamical processes act much faster than thermal or viscous ones, in the first approxima- 
tion one may consider only the dynamical structure of the disk. The disk dynamical structure is 
governed by three forces: gravity, which at a given place is fixed by the Kerr geometry; pressure; 
and rotation forces. The relative importance of each of these depend on the particular astrophysical 
situation — their values may either be "small" (dynamically unimportant) or "large" (dynamically 
important). Correspondingly, there are four "genuine" types of dynamical situations, as described 
in the upper half of Table [l] The four types of accretion disks we consider in this Living Review 
(thin, slim, thick, ADAF) are located in three of the four regions. The fourth region corresponds 
to a "free-fall" of dustj^j It is remarkable that the same types of accretion disks pop-up when 
one introduces other, very different, criteria for dividing the parameter space, such as accretion 
rate and opacity, as shown in lower half of Table [T] These four types of accretion disks are the 
next four Destinations on our road map. We characterize them in terms of: relative thickness, 
h = H/R; dimensionless accretion rate, m — O.lMc 2 / ZvEddj optical depth, t; importance of ad- 
vection, q — Q a dv/Qrad, where Q represents an energy flux; importance of radiation pressure, 
ft = -Pgas/ (-fgas + -Prad); location of inner edge, r; n ; and accretion efficiency rj. 

Destination 8: Thick Disk. 

h > 1, m > 1, r > 1, q ~ 1, < 1, r in ~ r mb , <q < 0.1 



8 This makes black hole accretion the most efficient energy generation process in the universe, short of matter- 
antimatter annihilation. 

9 Note that there is an analogy with stars, where there are three timescales - dynamical, thermal, and nuclear - 
which for most of the life of a star obey t^ yn <C t nuc . 

10 The slow rotation cases are sometimes referred to as "Bondi flows" in honor of Bondi's pioneering works on 
spherically-symmetric (non-rotating) accretion |50l . 
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Table 1: Upper half: The four a priori possible types of dynamical states of accretion structures, 
corresponding to the division of the dynamical parameter space into fast-slow rotation and large- 
small pressure. Lower half: A different division of the parameter space, corresponding to high-low 
accretion rate and large-small opacity. 





Fast rotation (Dis 


k) 


Slow rotation (Bondi) 


Large pressure 
Small pressure 


slim, thick 
thin 




ADAFs 
free-fall 




Accretion rate hi; 




Accretion rate low 


Large opacity 
Small opacity 


slim, thick 




thin 
ADAF 



Destination 9: Thin Disk. 

h < 1, m < 1, r > 1, q = 0, ~ 1, r in = r ms , rj ~ 0.1 

Destination 10: Slim Disk. 

h ~ 1, m > 1, t > 1, q ~ 1, < 1, r mb < r- m < r ms , rj < 0.1 

Destination 11: Advection- Dominated Accretion Flow (ADAF). 
h< 1, rh < 1, r < 1, q ~ 1, = 1, r mb < r in < r ms , i] < 0.1 

Interestingly, the parameter space of each of these types of accretion disks overlaps that of other 
solutions. For instance, ADAF solutions exist that have the same mass accretion rates as thin disk 
solutions 59 . In such cases, it is not clear how nature might choose one over the other. 

In Sections [|-[7j we will demonstrate that analytic (or at least semi-analytic) models exist for 
each of these types of accretion disks. However, this does not guarantee that they are stable. The 
issue of stability is our next destination. 

Destination 12: Stability. Stability analysis is important because the systematic differential ro- 
tation that is one of the defining characteristics of accretion disks is also a potential source 
of destabilizing energy. On the one hand, this may be essential, as the angular momentum 
transport and energy dissipation required for accretion may require disks to be mildly un- 
stable. On the other hand, if a model is violently unstable, then the basic assumption of a 
"steady-state" would be violated. 

Along with stability, we also look at the natural oscillation modes associated with accretion disks. 
The frequencies of some of these modes are tied to the properties of the black hole space-time, 
which is how these relate to the fundamental physics issues of interest in this Living Review. 

Destination 13: Oscillations. As with any finite distribution of fluid, accretion disks have nat- 
ural oscillation modes associated with them. If these modes can be excited at appreciable 
amplitudes, they may be able to modify the observed light curve of the disk in measur- 
able ways. This makes disk oscillations a leading candidate for explaining the quasi-periodic 



oscillations (QPOs) that we discuss in Section 12.4 



Because of the close observational links between black hole accretion disks and jets [92j|94], we 
include jets as the final Destination of our review. 

Destination 14: Jets. Jets are narrow (opening angle < 5°), long (length > 10 7 ly in the case of 
AGN), and fast (v > 0.9c) streams of matter emerging from very compact regions around the 
black hole, usually in opposite directions, presumably normal to the plane of the accretion 
disk. Jets can play a significant role in transporting energy and angular momentum away 



G 



from the accretion disk [JSj. They also play an important role in shaping the black hole's 
environment far beyond the gravitational reach of the black hole itself, affecting galactic 
evolution, particle acceleration, and intragalactic ionization. 

Going hand-in-glove with analytic models of accretion disks are direct numerical simulations. 
Although analytic theories have been extremely successful at explaining many general observational 
properties of quasars and microquasars, numerical simulations can be critically important in at least 
two respects: 1) as an extension of analytic work, by treating nonlinear perturbations and higher 
order coupling terms, and 2) in cases that are highly time variable or contain little symmetry, 
such that the prospects of finding an analytic solution are poor. There is also an important 
overlap region where various analytic and numerical methods are applicable and can be used to 
independently validate results. Because of these close connections between analytic and numerical 



work, we have a section (Section 111 dedicated to the discussion and review of direct numerical 



simulation of black hole accretion disks. 

We finish this Living Review with a section (Section 12 ) that tries to make some connections 
between the concepts discussed in earlier sections and actual observational phenomena. We em- 
phasize that we are not aiming to provide a comprehensive review of black hole observations, but 
rather to highlight a very small subset of these that are of particular relevance to our review. 

Throughout this review, we adopt the h ++ metric signature and often use units where 

c = 1 = G. To make all physical quantities dimensionless, we sometimes also use the mass of the 
black hole as a unit, M = 1. We use the common Einstein summation convention, where repeated 
indices in a formula imply summation over the range of that index. We also follow the common 
convention where Greek (Latin) indices are used for four-(three-)dimensional tensor quantities. 
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2 Three Destinations in Kerr's Strong Gravity 



In this section, we briefly describe the three destinations within Kerr's strong gravity that are most 
relevant to black hole accretion disk theory: 

1. Event Horizon: That radius inside of which escape from the black hole is not possible; 

2. Ergosphere: That radius inside of which negative energy states are possible (giving rise to 
the potentiality of tapping the energy of the black hole). 

3. Innermost Stable Circular Orbit (ISCO): That radius inside of which free circular 
orbital motion is not possible; 

Our principal question is: Could accretion disk theory unambiguously prove the existence of the 
event horizon, ergosphere, and ISCO using currently available or future observations? 

In realistic astrophysical situations involving astrophysical black holes (in particular quasars 
and microquasars) , the black hole itself is uncharged, and the gravity of accretion disk is practically 
negligible. This means that the spacetime metric g^ is given by the Kerr metric, determined by 
two parameters: total mass and total angular momentum J*. It is convenient to rescale them 
by 

M = ™ (2a) 
c z 

such that both M and a are measured in units of length. 

In the standard spherical Boyer-Lindquist coordinates the Kerr metric takes the form [3T] , 

/ 2Mr\ ff (r 2 + a 2 ) 2 -a 2 A sin 2 6 
= ^J' 9 = ?A ' 

2Mar sin 2 9 trb 2Mar 

9t<f> = 2 ' 9 



q 2 g 2 A ' 

2 2Ma 2 r sin 2 0\ . 2 . . A - a 2 sin 2 1 



9H>=\ T +« + ~2 S11T "> « - a 2 • 2, > 

Q J Ag z sin 

Q 2 2 rr ^ 68 1 fo\ 

X' 9ee = g , g g 3 

A g z g z 

2Mr + a 2 and g 2 = r 2 + a 2 cos 2 6. 
The Kerr metric depends neither on time t, nor on the azimuthal angle <j> around the sym- 
metry axis. These two symmetries can be expressed in a coordinate independent way by the two 
commuting Killing vectors rf = <5 p t and £ M = S^,, 

V^a + V A 77 M = 0, V^a+Va^ =0, (4a) 
W = fV/. (4b) 

Here V M denotes the covariant derivative, 

V^a = d^n x - T^r] v , = \ 2 g VK (dxg^n + d^g XK - d K g x ^) , (5) 

and <9 P = d/dx^ denotes the standard partial derivative. Formulae for the Kerr metric ([3| and all 
its non-zero Christoffell symbols |5j), are available from |305| . 
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In Boyer-Lindquist coordinates the t and <fi components of the Kerr metric can be expressed 
as scalar products of the Killing vectors, 

gtt = rfr)^ 9t<k = V 1 *^, 9U = fZi*- ( 6 ) 

From the Killing vectors one can also define the following constants of motion for a particle or 
photon with four-momentum p M 

energy: £ = -ifp^ = -p t , (7a) 
angular momentum: C = ^p^ = p$, (7b) 

specific angular momentum: I = — = — — = -, (7c) 

£ Pt ih 

Carter constant: C = K^p^Pu = (pef + t ■ (7d) 

sin 9 

The Carter constant C is connected to the Killing tensor K^„, which exists in the Kerr metric. 
Killing tensors obey, 

g [V«^ + VpKv* + V v K Rll + V K K V11 + V V K MK + V M Jf w ] - 0. (8) 

The other coordinates often used in black hole accretion disk research are the Kerr-Schild coordi- 
nates, in which the metric takes the form, 

9 M v = Vpu + fkufa, r) MV = diag(-l, 1, 1, 1) 
2Mr 3 

h 



ky 



r 4 
1. 


+ a 2 z 2 


rx 


+ ay 




+ a 2 ' 


ry 


— ax 




+ a 2 ' 


z 




r 





k z = -, (9) 
r 

where fc„ = (kt, k x ,k y ,k z ) is a unit vector, and r is given implicitly by the condition, 



4 ± 4+4=i- ( io ) 



Note that we have given the Kerr-Schild metric in its Cartesian form to prevent confusion with 
the spherical-polar Boyer-Lindquist coordinates. In keeping with this, unless specifically stated 
otherwise, the indices {t, 4>, r, 9} will always refer to the Boyer-Lindquist coordinates in this review. 



2.1 The event horizon 

The mathematically precise, general, definition of the event horizon involves topological consider- 
ations |207j . Here, we give a definition which is less general, but in the specific case of the Kerr 
geometry is fully equivalent. 

The Boyer-Lindquist coordinates split the Kerr spacetime into a "time" coordinate t and a 
three-dimensional "space," defined as t = const hypersurfaces. This split may be done in a coor- 
dinate independent way, based on the Killing vectors which exist in the Kerr spacetime. Indeed, 
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the family of non-geodesic observers A M with trajectories orthogonal to a family of 3-D spaces 
t = const is defined as, 



ZAMO: 
frame dragging: 



9t<t> 
94)4) 



(11a) 
(lib) 



They are called zero-angular-momentum-observers (ZAMO), because for them, the angular mo- 
mentum defined by ([7bl is zero, C = Na = 0. The ZAMO observers provide the standard of rest 



in the 3-D space: objects motionless with respect to the ZAMO frame of reference occupy fixed 
positions in space. 

We can also define a gravitational potential in the ZAMO frame: 



potential: 



$ = —In 
2 



= -2-1.9"! 



(12) 



The primary reason to call $ the gravitational potential is that, in Newton's theory, the observer 
who stays still in space experiences an acceleration due to "gravity" g^, which equals the gradient 
of the gravitational potential. In the Kerr spacetime it is, 



9 ii = M 



ZAMO 



N v V v Nu 



From (11a I one sees that at the surface 1/g = 0, the vector ry M is null, 



(13) 



0. Therefore, the 



ZAMO observers who provide the standard of rest, move on that surface with the speed of light. In 
order to stand still in this location, one must move radially out with the speed of lightj^] As it is 
clear from l/g u — is equivalent to A = 0. The last equation has a double solution, 



TG U±A 1 



a 

M 



(14) 



Note, that for r < r + the ZAMO "observers" are spacelike: standing still at a given radial location 
implies moving along a spacelike trajectory — i.e., faster than light. All trajectories that move 
radially out are also spacelike. Thus, the outer root r = r + of Eq. (14) defines the Kerr black 



hole event horizon: a null surface that surrounds a region from which nothing may escape. Outside 
the outer horizon (i.e., for r > r+) the normalization of A^ 1 is non-singular, and therefore the 



gravitational potential ( 12 ) is a non-singular, well-defined quantity. 



2.1.1 Detecting the event horizon 

One may think of two general classes of astrophysical observations that could provide evidence 
for a black hole horizon. Arguments in the first class are indirect; they are based on estimating a 
dimensionless "compactness parameter" 

(size of the object) 
(G/c 2 )(mass of the object) ' 

Arguments in the second class are more direct. They are based (in principle) on showing that some 
amount of radiation emitted by the source is lost inside the horizon. 

Evidence based on estimating the compactness parameter: A source for which observations 
indicate € ~ 1 may be suspected of having an event horizon. Values £ ss 1 have indeed been 

11 "Well, in our country, " said Alice, still panting a little, "you'd generally get to somewhere else if you run 
very fast for a long time, as we've been doing." "A slow sort of country!" said the Queen. "Now, here, you see, it 
takes all the running you can do, to keep in the same place. " 
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found in several astronomical sources. In order to know £, one must know mass and size of the 
source. The mass measurement is usually a direct one, because it may be based on an application 
of Kepler's laws. In a few cases the mass measurement is remarkably accurate. For example, in 
the case of Sgr A*, the supermassive black hole in the center of our Galaxy, the mass is measured 
to be M = (4.3 ± 0.5) x 10 6 M [TTI] . 

Until recently, estimates of size were always indirect, and generally not accurate. They are 
usually based on time variability or spectral considerations. For the former, the measurement rests 
on the logic that if the shortest observed variability time-scale is At, then the size of the source 
cannot be larger than R = cAt. For the latter, the argument goes like this: If the total radiative 
power L and the radiative flux F can be independently measured for a black-body source, then 
its size can be estimated from L — 4irR 2 F. Keep in mind that one must know the distance to the 
source in order to measure L. The flux can be estimated from F = aT A , where T is the temperature 
corresponding to the peak in the observed intensity versus frequency electromagnetic spectrum. 

It is hoped that in the near future, the next generation of high-tech radio telescopes will be able 
to measure directly the size of "the light circle" , which is uniquely related to the horizon size (see 
Figure [2]). For Sgr A*, at a distance of 8.28 ± 0.44 kpc |lllj . the event horizon corresponds to an 
angular size of ~ 10 /ias in the sky, making it an ideal target for near-future microarcsecond very 
long base interferometric techniques [77J IBS] • Here the plan is to measure the black hole "shadow" 
or "silhouette." 

Evidence based on the "no escape" argument: For accretion onto an object with a physical 
surface (such as a star), 100% of the gravitational binding energy released by accretion must be 
radiated away. This does not apply for a black hole since the event horizon allows for the energy 
to be advected into the hole without being radiated. This may allow for a black hole with an event 
horizon to be distinguished from another, similar-mass object with a surface, such as a neutron 
star. This argument was first developed by Narayan and collaborators 215|I216| 1214] : we describe 
it in more detail in Section 112.21 



2.2 The ergosphere 

In Newtonian gravity, angular momentum £ and angular velocity Q are related by the formula 
£ = r 2 fl, and therefore there is no ambiguity in defining a non-rotating frame as Q = = I. 
However, in the Kerr geometry £ oc (CI — u>), where uj = —gt<f>/ 9if>4> is the angular velocity of the 
frame dragging induced by the Lense-Thirring effect. Therefore, O = does not imply £ = 0. This 
leads to two different standards of "rotational rest": the Zero Angular Velocity Observer (ZAVO) 
and the Zero Angular Momentum Observer (ZAMO), 

ZAVO frame {fl = 0): = {-rfr) v )- 1/2 V 1 *, (16) 

ZAMO frame (£ = 0): N» = e*^ + u^) . (17) 

These two frames rotate with respect to each other with the frame-dragging angular velocity 

w = —9t<f>l '94><j>- 

The ZAMO frame defines a local standard of rest with respect to the local compass of inertia: a 
gyroscope stationary in the ZAMO frame does not precess. Considering its kinematic invariants J^j 

12 For a congruence of observers (or particles or photons) with four velocity f/ M , the kinematic invariants fully 
describe their relative motion. Consider those that, in a particular moment so, occupy the surface of an infinitesi- 
mally small sphere. Now, consider the deformation of that surface at a later moment so + ds. The volume change 
dV/ds = is called expansion. The shear tensor u^ Ll/ measures the ellipsoidal distortion of this sphere, and the 
vorticity tensor w M „ describes its rotation (i.e., three independent rotations around three perpendicular axes). Ex- 
pansion, shear, and vorticity are determined by the tensor X M „ = V ^U v in the following way: = (1/3)JC„, 
tr^ = h a ^ bP v (l/2)(X aj3 + Xp a ) -eh» u , and = h\ h f i ({l/2)(X af} - Xp a ). Here /»« = <5 Q M + U a U„ is the 
projection tensor. The acceleration a^ is also considered a kinematic invariant. 
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Figure 2: Silhouettes of Sgr A* calculated for four optically thin accretion structures, characterized 
by very different physical conditions. The display is intentionally reversed in black-and-white and 
saturated in order to better show the less luminous parts. Although "dirty astrophysics" makes 
the most prominent differences, effects of the "pure strong gravity" are also seen in the form of 
"the light circle" , a tiny almost circular feature at the center. Its shape and size depends only on 
the black hole mass and spin. Image reproduced by permission from [297], copyright by ESO. 
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one sees that the ZAMO frame is non-inertial (a„ 7^ 0), non-rigid (a^ 7^ 0, = 0), and surface- 
forming (u)^ = 0). The ZAMO vectors ff and _/V M are time-like everywhere outside the horizon, 
i.e., outside the surface l/g tt = 0. This means that the energy of a particle or photon with a 
four-momentum p M measured by the ZAMO is positive, £?zamo = N^p^ > 0. 

The ZAVO frame defines a global standard of rest with respect to distant stars: a telescope that 
points to a fixed star does not rotate in the ZAVO frame. Considering its kinematic invariants one 
sees that the ZAVO frame is non-inertial (a„ 7^ 0), rigid {(Jav = 0, = 0), and not surface- forming 
(utfiu 7^ 0). At infinity, i.e., for r — > 00, it is (r] v ri^) — g tt — > —1, and therefore rt L — > For this 
reason, rf L is called the stationary observer at infinity. The ZAVO vectors rf 1 and are timelike 
outside the region surrounded by the surface gu — 0, called the ergosphere. Inside the ergosphere 
77^ and are spacelike. This means that inside the ergosphere, the conserved energy of a particle 



(i.e., the energy measured "at infinity"), as defined by (7a), may be negative. 

Penrose |242j considered a process in which, inside the ergosphere, a particle with energy 
E°° > decays into two particles with energies > and E°° = — \E°°\ < 0. The particle with 
positive energy escapes to infinity, and the particle with the negative energy gets absorbed by the 
black hole. Then, because = E — E°° = E + \E°°\ > E, one gets a net gain of positive energy 
at infinity. The source of energy in this Penrose process is the rotational energy of the black hole. 
Indeed, the angular momentum absorbed by the black hole, J°° = pig* is necessarily negative (in 
the sense that J°°ujh < 0), which follows from 

£zamo = -ft e*(ry* + ) = e*(£~ - uj h J°°) > 0, 

and thus 

lu h J°° < E x < 0. (18) 

A more complete presentation of the Penrose process is made in |310j . At this time it appears the 
most likely realization of the Penrose process would be the Blandford-Znajek mechanism [JS] for 
launching jets from quasars and microquasars. Observations suggest [255, 252 , and simulations 
confirm 303, 304 , that through this mechanism it is possible to extract more energy from the 
system than is being delivered by accretion. We discuss jets and the Blandford-Znajek mechanism 
more in Sections [TOl and fll.71 

2.3 ISCO: the orbit of marginal stability 

Particles (with velocity normalization u^u^ = —1) and photons (with velocity normalization 
u^u^ — 0) move freely on "geodesic" trajectories x M = x M (s), with velocities it M = dx^/ds, 
characterized by vanishing accelerations 

d 2 r u dr K dr^ 

A constant of motion X, such as those defined in Eq. is conserved along a geodesic trajectory 



(19) in the sense that u^V = 0. 

Circular geodesic motion in the equatorial plane (9 — 7r/2) is of fundamental importance in 
black hole accretion disk theory. The four velocity corresponding to circular motion is defined by, 

u" =A(rj IJ -+n^) ) (20) 

where f2 = u^/u 1 = dcf>/dt is the angular velocity measured by the stationary observer (ZAVO, see 



Section 2.2), and the redshift factor, A = u*, follows from u^u v g^ u = — 1, 

-A- 2 = g tt + 2n gt4> + n 2 g4>4> . (21) 
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Other connections between these quantities that are particularly useful in our later calculations 
also follow from u^u v g^ u = — 1: 

n (■ 9tt + gt<p e n g^ + gt<f> A , n v , 00 s 

It is convenient to define the effective potential, 

W off = -iln| ff "-2V* + ^ 2 5*1, (23) 

because in terms of W e ff and the rescaled energy £* = ln£, slightly non-circular motion, i.e., with 
V 2 = u r u r g rr + u°u e ggg <C u^u^g^, is characterized by the equation, 

^V 2 = S*-U eS , (24) 

which has the same form and the same physical meaning as the corresponding Newtonian equation. 
Therefore, exactly as in Newtonian theory, unperturbed circular Keplerian orbits are given by the 
condition of an extremum (minimum or maximum) of the effective potential (6 = tt/2), 

.(^)r°W(£)-*(£) +<, (£)-°}- 

This quadratic equation for I has two roots I — ±^j<-(r, a), corresponding to "corotating" and 



"countcrrotating" Keplerian orbits. Their explicit algebraic form is given in Eq. (35) in Section 2.5 
As in Newtonian theory, slightly non-circular orbits (with V ^ being either Sr or 89) are fully 
determined by the simple harmonic oscillator equations, 

Sf + w, 2 Sr = 0, 59 + u 2 59 = 0, (26) 

where the radial w r and vertical u>g epicyclic frequencies are second derivatives of the effective 
potential, 



where dxj 2 = —g xx dx 2 . The epicyclic frequencies (27 1 are measured by the comoving observer. 
To get the frequencies w* measured by the stationary "observer at infinity" (Section |2.2| , one must 
rescale by the redshift factor ui* — Auj x . Obviously, when (u> r ) 2 < 0, the epicyclic radial oscillations 



described by Eq. ( 26 ) are unstable - from Eq. ( 27 1 we see that they correspond to maxima of the 



effective potential. This happens for all circular orbits with radii less than r = r ms (a), and this 
limiting radius is called ISCO, the innermost stable circular orbit. 

Free circular orbits with r > r ms are stable, while those with r < r ms are not. Accordingly, 
accretion flows of almost free matter (i.e., with stresses insignificant in comparison with gravity 
or centrifugal effects), resemble almost circular motion for r > r ms , and almost radial free-fall for 
r < r ms- For thin disks, this transition in the character of the flow is expected to produce an 



effective inner truncation radius in the disk (see Section 5.3). The exceptional stability of the inner 
radius of the X-ray binary LMC X-3 |293j . provides considerable evidence for such a connection 
and, hence, for the existence of the ISCO. The transition of the flow at the ISCO may also show 
up in the observed variability pattern, if variability is modulated by the orbital motion. In this 
case, one may expect that the there will be no variability observed with frequencies v > z^iscOj 
i.e., higher than the Keplerian orbital frequency at ISCO, or that the quality factor for variability, 
Q ~ v I ISv will significantly drop at i^isco- Several variants of this idea have been discussed [331134) . 
and some observational evidence to support them has been presented (see Figure [3]). 
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Figure 3: Evidence for the existence of the ISCO from data recorded by the Rossi X-ray Timing 
Explorer satellite from neutron star binary source 4U 1636-536 |33J. The source shows quasi- 
periodic oscillations (QPOs) with frequencies in the range 650 Hz < v < 900 Hz. The sharp drop 
in the quality factor (bottom panel) seen at ~ 870 Hz may be attributable to the ISCO [34] . 



2.4 The Paczyriski— Wiita potential 

For a non-rotating black hole (a = 0), the Kerr metric reduces to the Schwarzschild solution, 

_ d t 2 + U- 2 J£^j 1 dr 2 + r 2 [d9 2 + sin 2 9 d(f> 2 ] . (28) 

Paczyhski and Wiita [236] proposed a practical and accurate Newtonian model for a Schwarzschild 
black hole, based on the gravitational potential, 

GM 2GM 
®pw = , r s = ^^. 29 

r — rg &■ 

The Paczynski-Wiita potential became a very handy tool for studying black hole astrophysics. It 
has been used in many papers on the subject and still has applicability today. The Schwarzschild 
and the Paczyhski- Wiita expressions for the Keplerian angular momentum and locations of the 



marginally stable and marginally bound orbits (Section 2.3) are identical. Similar, though less 
commonly adopted, pseudo-Newtonian potentials have also been found for Kerr (rotating) black 
holes |251 IZ7S1I2I3] . 



2.5 Summary: characteristic radii and frequencies 

We end this section with a few formulae for the Kerr geometry that we will use elsewhere in this 
review. 

Keplerian circular orbits exist in the region r > r p h, with r p h being the circular photon orbit. 
Bound orbits exist in the region r > r m b, with r m b being the marginally bound orbit, and stable 
orbits exist for r > r ms , with r ms being the marginally stable orbit (also called the ISCO - 
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Section 2.3). The location of these radii, as well as the location of the horizon ru and ergosphere 
ro, are given by the following formulae |31j : 



photon r p h 

bound r m b 
stable r ms 
horizon rjj 
ergosphere rg 



2r r , < 1 + cos 



* -if \ 
-cos (a*) 



2r G (l - y + \/l-a* 
r G {3 + Z 2 - [(3 - Z x )(3 + Zi + 2Z 2 )] 1/2 } 

r G (l + v/l - a 2 cos 2 , 



(30) 

(31) 
(32) 
(33) 
(34) 



where Z x = l+tl-a^/^l+a^+Cl-a,) 1 / 3 ], Z2 = (3a 2 + Z 2 ) i/2 ; ^ = fl/Mj and Tq = GM/c 2 
is the gravitational radius. 

The Keplerian angular momentum ix and angular velocity 91 k, and the angular velocity of 
frame dragging u are given by, 



-K 



O 2 



M l ' 2 (r 2 -2aM 1 / 2 r 1 / 2 + a 2 ) 
' K ' ' ' " ' ~ r 3/2 _ 2MrV2 + aMV 2 



GM 



(r^+aM 1 / 2 ) 
2aMr 



(a 2 +r 2 )g 2 + 2a 2 Mr sin 2 9 ' 
The epicyclic frequencies measured "at infinity" are (here x = r/M), 

K) 2 = (}% (l - 6X' 1 + 8a*x~V 2 - 3a 2 x- 2 ) 



*\2 



K) 



9 2 



2 K (l-4a*x~ 3/2 +3alx~ 2 ^j 



(35) 

(36) 

(37) 
(38) 

(39) 
(40) 



Comparing the Keplerian and epicyclic frequencies and the characteristic radii between the Schwarz- 
schild metric and the Paczyhski-Wiita potential (Section 2.4), we find for the Schwarzschild metric, 



9. 



GM 



K 



K) 2 , (w^n^i-ear 1 ), 



and x ms = 6, x mh = 4, x ph = 3, 
and for the Paczyhski-Wiita potential, 



(41) 



O 2 - 

\L K — 



GM 



(1 - 62;- 1 ) 



3(l-2.x-i) 2 "" K) ' K * } - ^(1-21-1)' 



and x ms = 6, x mh = 4. 



(42) 
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3 Matter Description: General Principles 



Having provided a detailed description of the key signatures associated with a black hole spacetime, 
we now move into the mirkier realm of the accretion disk itself. We start from the fundamental 
conservation laws that govern the behavior of all matter, namely the conservation of rest mass and 
conservation of energy-momentum, stated mathematically as 

V„(pO = 0, V^ = 0. (43) 

Here p is the rest mass density, is the four velocity of matter, and T£ is the stress energy 



tensor describing properties of the matter. The conservation equations (43) are supplemented 
by numerous "material" equations, like the equation of state, prescriptions of viscosity, opacity, 
conductivity, etc. Several of them arc phcnomcnological or simple approximations. Nevertheless, 
we can give a GEN-eral form of T£ that is relevant to accretion disk theory as a sum of FLU-id, 
VIS-cous, MAX-well, and RAD-iation parts, which may be written as, 

rad, (44) 
(45) 
(46) 



(T^)gen 


= CC)flu + 


CC)vis + (2?)max + CC) 


CC)flu 


= (pu^iWu 


v ) + snp, 


W)vis 


= v*<r2, 


i^F^F^, 




= F» a F av - 


(^T)rad 


4 


-u^Fv+Uv F". 



(47) 
(48) 

Here W = enthalpy, 6£ — Kronecker delta tensor, P = pressure, v* = kinematic viscosity, 
= shear, F^ v = Faraday electromagnetic field tensor, E = radiation energy density, and F' 1 = ra- 
diation flux. In the remainder of this section we describe these components one by one, including 
the most relevant details. Most models of accretion disks are given by steady-state solutions of the 



conservation equations (43), with particular choices of the form of the stress-energy tensor Tjf, and 
a corresponding choice of the supplementary material equations. For example, thick accretion disk 
models (Section [4| often assume (T^)vis = (T^)max = (F^)rad = 0, thin disk models (Section [5]) 



assume (T^')max = 0, and most current numerical models (Section 11) assume (T^)rad = 



3.1 The fluid part 

The one absolutely essential piece of the stress-energy tensor for describing accretion disks is the 
fluid part, (X^)flu = {pu tl ){Wu v ) + 5% P. The fluid density, enthalpy, and pressure, as well as 
other fluid characteristics, are linked by the first law of thermodynamics, dll — T dS — P dV , which 
we write in the form, 

de = Wdp + nTdS, (49) 

where U is the internal energy, T is the temperature, S is the entropy, and e = pc 2 + H is the total 
energy density, with n being the internal energy density, and 

V=K 17=5, W =Z±±. (50 ) 
n n p 

The equation of state is often assumed to be that of an ideal gas, 

P = -pT, (51) 
with TZ being the gas constant and p: the mean molecular weight. 
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Sometimes we may wish to consider a two temperature fluid, where the temperature T t and 
molecular weight //, of the ions are different from those of the electrons (T e and /x e ). For such a 
case 

P = Pi + P e = —pT, + —pT e . (52) 

P~i Me 

Two-temperature plasmas are critical in advection-dominated flows (discussed in Section [7]). Two- 
temperature fluids are also important when one considers radiation [281] . as it is the ions that are 
generally heated by dissipative processes in the disk, while it is generally the electrons that radiate. 
Ions and electrons normally exchange energy via Coulomb collisions. As this process is generally 
not very efficient, the electrons in the inner parts of accretion flows are usually much cooler than 
the ions (Coulomb collisions are not able to heat the electrons as fast as they radiate or advect into 
the black hole). However, there have been suggestions that more efficient processes may couple 
the ions and electrons |245j . such as plasma waves [35] or kinetic instabilities [282) . At this point, 
this remains an open issue in plasma physics, so it is difficult to know how much heating electrons 
experience. 



3.1.1 Perfect fluid 



In the case of a perfect fluid, the whole stress-energy tensor (|44j) is given by its fluid part (45), and 
all other parts vanish, i.e., (T^)gen = (T^)flu- In this particular case, one can use V ^(T^t]") = 
r/ u V^(T>f) + T^V^) = + = 0, and similarly derived V^T^f) = 0, to prove that 



B = -W{uri) = -W u u J = W(u£) = W U4> , 



(53) 



are constants of motion. We can identify B as the Bernoulli function and J as the angular 
momentum. Their ratio is obviously also a constant of motion, 



B 



(54) 



identical in form with the specific angular momentum (7c), which is a constant of geodesic motion. 



3.2 The stress part 



In the stress part (T^)vis = v * the shear tensor is a kinematic invariant (cf. Footnote 12) 
It is defined as 

"1 



' jJLU 



(55) 



where the symbol j_ denotes projection into the instantaneous 3-space perpendicular to u M in the 
sense that {X^)± = X a (8^ a + u a ). The other kinematic invariants are vorticity, 



<>nv 



(56) 



and expansion, 



(57) 



In the standard hydrodynamical description (e.g. |168| ). the viscous stress tensor, Sff, is pro- 
portional to the shear tensor, 

S£ = v.p*Z. (58) 
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The rate of heat generation by viscous stress in a volume V is then given by 

Q+ = J S»a\dV. (59) 

In addition, the rates of viscous angular momentum and energy transport across a surface S, with 
a unit normal vector iV, are 

Js= I ' S2S k v m NndS, B s = j SWN^dS. (60) 

For the case of purely circular motion, where u l = A{rj l + f2£ J ), the kinematic invariants are 

e = o, o-^ = iA 8 * 2 a M n, u^ = ^A 3 ^-%t, (6i) 

where ^ 2 = g\± — gttg<f>4>- It is a general property that (X^)±u^ — 0, and so for purely circular 
motion, one has, 

o^vrf = -ficr^", u^*/" = '. (62) 



From Eqs. (60) and (62) one deduces that for purely circular motion, the rates of energy and 



angular momentum transport are related as 

B s = -n s J s , (63) 

where fig is the angular velocity averaged on the surface S. From this, one sees that as angular 
momentum is transported outward, additional energy is carried inward by the fluid. 

3.2.1 The alpha viscosity prescription 

As we mentioned in Section [T] the viscosity in astrophysical accretion disks can not come from 
ordinary molecular viscosity, as this is orders of magnitude too weak to explain observed phenom- 
ena. Instead, the source of stresses in the disk is likely turbulence driven by the magneto-rotational 



instability (MM, described in Section 8.2). Even so, one can still parametrize the stresses within 
the disk as an effective viscosity and use the normal machinery of standard hydrodynamics with- 
out the complication of magnetohydrodynamics (MHD). This is sometimes desirable as analytic 
treatments of MHD can be very difficult to work with and full numerical treatments can be costly. 



For these reasons, the Shakura-Sunyaev "alpha viscosity" prescription [279 still finds applica- 
tion today. It is an ad hoc assumption based on dimensional arguments. Shakura and Sunyaev 
realized that if the source of viscosity in accretion disks is turbulence, then the kinematic viscosity 
coefficient has the form, 

v* ^l v 0l (64) 

where lo is the correlation length of turbulence and vq is the mean turbulent speed. Assuming that 
the velocity of turbulent elements cannot exceed the sound speed, vq < cs, and that their typical 
size cannot be greater than the disk thickness, Iq < H, one gets 

= a Hcs , (65) 

where < a < 1 is a dimensionless coefficient, assumed by Shakura and Sunyaev to be a constant. 
For thin accretion disks (see Section u5|) the viscous stress tensor reduces to an internal torque 



with the following approximate form [see Eqs. (55) and (58)] 



Trci, « pv*r— . (66) 
or 
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However, for thin disks, r(dQ/dr) w — Q and c$ ~ (P/p) 1 / 2 ~ SIP, so Shakura and Sunyaev argued 
that the torque must have the form T r <j> — —aP. A critical question that was left unanswered was 
what pressure P one should consider: P ga s, P-ad, or Plot = Pgas + Pad? This question has now 
been answered using numerical simulations [128], so that we now know the appropriate pressure 
to be Pxot- Typical values of a estimated from magnetohydrodynamic simulations are close to 
0.02 [122] . while observations suggest a value closer to 0.1 (see jl48j and references therein). 



3.3 The Maxwell part 

Magnetic fields may play many interesting roles in black hole accretion disks. Large scale magnetic 
fields threading a disk may exert a torque, thereby extracting angular momentum [48] . Similarly, 
large scale poloidal magnetic fields threading the inner disk, ergosphere, or black hole, have been 
shown to be able to carry energy and angular momentum away from the system, and power jets |49) . 
Weak magnetic fields can tap the differential rotation of the disk itself to amplify and trigger an 
instability that leads to turbulence, angular momentum transport, and energy dissipation (exactly 
the processes that are needed for accretion to happen) (251 HZ] • 

In most black hole accretion disks, it is reasonable to assume ideal MHD, whereby the conduc- 
tivity is infinite, and consequently the magnetic diffusivity is zero. Whenever this is true, magnetic 
field lines are effectively frozen into the fluid. A corollary to this is that parcels of fluid are re- 
stricted to moving along field lines, like "beads" on a wire. In ideal MHD, the Faraday tensor 
obeys the homogeneous Maxwell's equation 

V„(*F^) = 0, (67) 

where *Fjf is the dual. If we define a magnetic field 4- vector b^ = u v *F^, then using b^u^ = 
one can show that 

*F>t = Vu v - b v u^. (68) 



Using this, it is easy to show that the spatial components of (67) give the induction equation 

d t {^-gB l ) = -d,[V^(BV - B j v% (69) 
while the time component gives the divergence-free constraint 

d i ( y /=gB i ) = o, (70) 

where B l — u l b l — u 1 ^, and g is the 4-metric determinant. 



3.3.1 The magneto-rotational instability (MRI) 



We mentioned in Section 3.2 that a hydrodynamic treatment of accretion requires an internal 
viscous stress tensor of the form < 0. However, we also pointed out that ordinary molecular 
viscosity is too weak to provide the necessary level of stress. Another possible source is turbulence. 
The mean stress from turbulence always has the property that T r ^ < 0, and so it can act as an ef- 
fective viscosity. As we will explain in Section [8.2[ weak magnetic fields inside a disk are able to tap 
the shear energy of its differential rotation to power turbulent fluctuations. This happens through 
a mechanism known as the magneto-rotational (or "Balbus-Hawley" ) instability II 181 ITT] . Al- 
though the non-linear behavior of the MRI and the turbulence it generates is quite complicated, its 
net effect on the accretion disk can, in principle, be characterized as an effective viscosity, possibly 
making the treatment much simpler. However, no such complete treatment has been developed at 
this time. 
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3.4 The radiation part 

Radiation is important in accretion disks as a way to carry excess energy away from the system. 



In geometrically thin, optically thick (Shakura-Sunyaev) accretion disks (Section 5.3), radiation is 
highly efficient and nearly all of the heat generated within the disk is radiated locally. Thus, the 
disk remains relatively cold. In other cases, such as ADAFs (Section [7]), radiation is inefficient; 
such disks often remain geometrically thick and optically thin. 

In the optically thin limit, the radiation emissivity / has the following components: brems- 
strahlung /b r , synchrotron / sy nch, and their Comptonized parts /b r ,c an d /synch, c- In the optically 
thick limit, one often uses the diffusion approximation with the total optical depth r = r a b s + t cs 
coming from the absorption and electron scattering optical depths. In the two limits, the emissivity 
is then 

, ( /br + /synch + ,/br.c + /synch,c optically thin case (r < 1), 

J = 1 8erT 4 \' / 

I 3H ^ optically thick case (t> 1), 

where a is the Stefan-Boltzmann constant. In the intermediate case one should solve the transfer 
equation to get reliable results, as has been done in [288, ]5B]. Often, though, the solution of the 
grey problem obtained by Hubeny |133) can serve reliably: 



/ = 



AaT, 



4 



— + H (/br + /synch + /br,C + /synch, c) 1 



(72) 



In sophisticated software packages like BHSPEC, color temperature corrections in the optically thick 
case (the "hardening factor") are often applied |66| . 

In the remaining parts of this section we give explicit formulae for the bremsstrahlung and syn- 
chrotron emissivities and their Compton enhancements. These sections are taken almost directly 
from the work of Narayan and Yi [225j . Additional derivations and discussions of these equations 
in the black hole accretion disk context may be found in [599 , 295 , 225 , [57] . 

3.4.1 Bremsstrahlung 

Thermal bremsstrahlung (or free-free emission) is caused by the inelastic scattering of relativistic 
thermal electrons off (nonrelativistic) ions and other electrons. The emissivity (emission rate per 
unit volume) is /b r = fei + fee- The ion-electron part is given by [225] 

2 f4(2%) 1/2 (l + 1.7810 e 134 ) 9 e <l, 
f ei ; = n e na T ca f m e c 2 x < n } * 3 > K e ' e ' (73) 

' \f^ln(1.1230 e +0.48) + 1.5] e > 1. 

where n e is the electron number density, n is the ion number density averaged over all species, 
<tt = 6.62 x 10 - 25 cm 2 is the Thomson cross section, a/ = 1/137 is the fine structure constant, 
S e = kBT e /m e c 2 is the dimensionless electron temperature, and ks is the Boltzmann constant. 
The electron-electron part is given by |225j 

2 2 2 1 ^(44 - 37r 2 )^ /2 (l + l.l(9 e + 2 - 1.250c 5/2 ) 9 e < 1, 
f = nicrim e c af x < 97r ' y ' K e ' (74) 

J e e ; [24f9 e [ln(1.1232(9 e ) + 1.28] 9 e > 1. 

where r e = e 2 /m e c 2 is the classical radius of the electron. 

3.4.2 Synchrotron 

Assuming the accretion environment is threaded by magnetic fields, the hot (relativistic) electrons 
can also radiate via synchrotron emission. For a relativistic Maxwellian distribution of electrons, 
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the formula is [225] 



/synch" 3c 2 K1 - dr 



47rm„c 



(75) 



where e is the electric charge, B is the equipartition magnetic field strength, and xm is the solution 
of the transcendental equation 

ex P (1.88994f) = 2.49 x l<r"» (^f) ^ + ^ f^^^ , (*» 

where the radius r must be in physical units and K 2 is the modified Bessel function of the 2nd 
kind. This expression is valid only for 9 e > 1, but that is sufficient in most applications. 

3.4.3 Comptonization 

The hot, relativistic electrons can also Compton up-scatter the photons emitted via bremsstrahlung 
and synchrotron radiation. The formulae for this are |225j 

, , / mxc % [3-("3+D - ( 3 g e )-fa+i)] | 

/synch,C = /synch [r]l ~ V2 (x c /9 e ) n3 } . (78) 

Here 77 = 1 + 771 + 772(2; / '8 e)'' 3 is the Compton energy enhancement factor, and 

hv hv c Xi(x\ — 1) 



m e c z m e c z 1 — X1X2 

Xl = i + 4e e + we 2 e , "n = ^, 

X 2 = 1 - exp (—Tea), r] a = —1 — ]xix 2 / Inxi, (79) 

where h is Planck's constant and v c is the critical frequency, below which it is assumed that the 
emission is completely self-absorbed and above which the emission is assumed to be optically thin. 
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4 Thick Disks, Polish Doughnuts, &; Magnetized Tori 



In this Section we discuss the simplest analytic model of a black hole accretion disk - the "Polish 
doughnut." It is simplest in the sense that it only considers gravity (Section[2]), plus a perfect fluid 



(Section |3.1.1 1, i.e., the absolute minimal description of accretion. We include magnetized tori 
in Section |4.2 which allows for (T^) M ax ^ 0, but otherwise (T^) V is = (T^)max = (^)rad = 
throughout this section. 

4.1 Polish doughnuts 

Paczyhski and his collaborators developed, in the late 1970s and early 1980s, a very general method 
of constructing perfect fluid equilibria of matter orbiting around a Kerr black hole |139l 12361 12351 
1234] . They assumed for the stress energy tensor and four velocity, 

n = (Tt) FLU = pwu^+ps^ 

u" - A(rf + n^), (80) 



and derived from V M T^ = that, 



V v \nA-^= l V v P. (81) 
1 — til p 



In the case of a barytropic fluid P = P(e), the right hand side of Eq. (81) is the gradient of a 
scalar function, and thus the left hand side must also be the gradient of a scalar, which is possible 
if and only if 

l = i(Sl). (82) 

This statement is one of several useful integrability conditions, collectively called von Zeipel theo- 
rems, found by a number of authors [5TJ (5^1 [TJ I156| . 

In real flows, the function £ = £(il) is determined by dissipative processes that have timescales 
much longer than the dynamical timescale, and are not yet fully understood. Paczyhski realized 
that instead of deriving I — £(Q) from unsure assumptions about viscosity that involve a free 
function fixed ad hoc (e.g., by assuming a(r,6) = const), one may instead assume the result, i.e., 
assume i — £(Q). Assuming £ = £(Q) is not self-consistent, but neither is assuming a(r, 9) = const. 

In Boyer-Lindquist coordinates, the equation for the equipressure surfaces, P — P(r, 9) = 
const, may be written as rp = rp(6), with the function rp(6) given by 

drp = 8e_P = {l~£n)d e \nA + £d e n 

d6 d r P (l-£n)d r \n A + £d r n' [ ' 

Using the expressions for A = A(r, 6), fl — Q(r,9), and £ = £{r,9) from Section [2] (Eqs. 21 and 
22), one can integrate Eq. (83) to get the equipressure surfaces. A description of how to do this 
for both Schwarzschild and Kerr black holes is given in [57]. Figure [4] illustrates the simplest (and 
important) case of I = £(Sl) = £q = const. 

Another useful way to think about thick disks is from the relativistic analog of the Newtonian 
effective potential $, 

r p dP 

<&-$!„ = -/ — , (84) 



/o pW 

where $; n is the potential at the boundary of the thick disk. For constant angular momentum £, 
the form of the potential reduces to $ = ln(— Ut). Provided £ > £ ms , the potential $(r, 9) will have 
a saddle point $ C us P at r — r cusp , 9 = n/2. We can define the parameter A$ = <t> in — $ cusp as 
the potential barrier (energy gap) at the inner edge of the disk. If A$ < 0, the disk lies entirely 
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Figure 4: In equilibrium, the equipressure surfaces should coincide with the surfaces shown by 
the solid lines in the right panel. Note the Roche lobe, self-crossing at the cusp. The cusp 
and the center, both located at the equatorial plane 9 = n/2, are circles on which the pressure 
gradient vanishes. Thus, the (constant) angular momentum of matter equals the Keplerian angular 
momentum at these two circles, £q — ^/f ( r cusp) = &k (/center), as shown in the upper left panel. 
Image adapted from [98]. In this figure W refers to the effective potential. 



within its Roche lobe, whereas if A<& > 0, matter will spill into the black hole even without any 
loss of angular momentum. 

Before leaving the topic of Polish doughnuts, we should point out that, starting with the work 
of Hawley, Smarr, and Wilson |125] . this simple, analytic solution has been the most commonly 
used starting condition for numerical studies of black hole accretion. 



4.2 Magnetized Tori 

Komissarov |156j was able to extend the Polish doughnut solution by adding a purely azimuthal 
magnetic field to create a magnetized torus. This is possible because a magnetic field of this form 
only enters the equilibrium solution as an additional pressure-like term. For example, the extended 



form of Eq. ( |8T| ) is 

iv„n = YjdL + v ^ 2p m a g ) 
i- in p * 2 p 



V, lnA-f^ = ^ + v ^*?"«' , (85) 



where V 2 = g t <t>gt4> - gttg</></> and Eq. (|84j) becomes 

(86) 

where P mag = ^Prnag- Komissarov [156] gives a procedure for solving the case of a barotropic 
magnetized torus with constant angular momentum (£ — const.). 
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5 Thin Disks 



Most analytic accretion disk models assume a stationary and axially symmetric state of the matter 
being accreted into the black hole. In such models, all physical quantities depend only on the two 
spatial coordinates: the "radial" distance from the center r, and the "vertical" distance from the 
equatorial symmetry plane z. In addition, the most often studied models assume that the disk is 
not vertically thick. In "thin" disks z/r <C 1 everywhere inside the matter distribution, and in 
"slim" disks (Section [6]) z/r < X. 

In thin and slim disk models, one often uses a vertically integrated form for many physical 
quantities. For example, instead of density p(r, z) one uses the surface density defined as, 

r+H(r) 

S(r) = / p(r,z)dz, (87) 
where z = ± H(r) gives the location of the surface of the accretion disk. 



5.1 Equations in the Kerr geometry 

The general relativistic equations describing the physics of thin disks have been derived indepen- 
dently by several authors |169[ [7j [T3l H06| [40] . Here we present them in the form used in [268J : 
(i) Moss conservation (continuity): 

M = -2ttSA 1/2 - V 

Vl - V 2 

where V is the gas radial velocity measured by an observer at fixed r who co-rotates with the fluid, 
and A has the same meaning as in Section [2] 



(ii) Radial momentum conservation: 

V dV A 1 dP 



1 - V 2 dr r £ dr 



(89) 



where 



A= ma 

r 3 An+n K i-n 2 R 2 ' 1 ; 

A — (r 2 + a 2 ) 2 — a 2 Asin 2 6*, f2 = u^/u 1 is the angular velocity with respect to the station- 
ary observer, Cl = fl — to is the angular velocity with respect to the inertial observer, £1^ = 
±M 1 / 2 /(r 3 / 2 ± aM 1 / 2 ) are the angular frequencies of the co-rotating and counter-rotating Keple- 
rian orbits, and R — A/(r 2 A 1 / 2 ) is the radius of gyration, 
(iii) Angular momentum conservation: 

M i 1 / 2 A 1 / 2 7 

~r-A«) = -«n (91) 

Zir r 

where C = is the specific angular momentum, 7 is the Lorentz factor, II = 2HP can be 



considered to be the vertically integrated pressure, a is the standard alpha viscosity (Section 3.2.1 1, 
and Ci n is the specific angular momentum at the horizon, which can not be known a priori. As we 
explain in the next section, it provides an eigenvalue linked to the unique eigensolution of the set 
of thin disk differential equations, once they are properly constrained by boundary and regularity 
conditions. 

(iv) Vertical equilibrium: 

n _ c?-a\e*-i) 
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with £ = —Ut being the conserved energy associated with the time symmetry, 
(v) Energy conservation: 

aUAj 2 dQ 32ctT 4 _ M 1 ( dP r ^dp\ 



r 3 dr 3 kS 2irrp T3 — 1 \ dr p dr J 

where T is the temperature in the equatorial plane, k is the mean (frequency-independent) opacity, 

Ti =0* + (4-3/3*)(r 3 -l), 

r =1 , (4-3^)( 7g -l) 
3 + 12(l-/3//3 m )( 7ff -l)+/3' 

/3 = Pgas/(-P g as + Prad + Pmag), Pm = fgas/C-Pgas + ^mag), /?* = /3(4 - /3 m )/3/3 m , and 7g is the ratio 
of specific heats of the gas. 

5.2 The eigenvalue problem 

Through a series of algebraic manipulations one can reduce the thin disk equations to a set of two 
ordinary differential equations for two dependent variables, e.g., the Mach number Ai = —Vjcg — 
—VTi/P and the angular momentum C = ua,. Their structure reveals an important point here, 

dlnTW _ A/~i(r,A4,£) dln£ _ Af 2 (r,M,C) 
dlnr ~ V(r,M,C) dlnr ~ V(r,M,£) ' 

In order for this to yield a non-singular physical solution, the numerators Afi and A2 must vanish 
at the same radius as the denominator T>. The denominator vanishes at the "sonic" radius r son j C 
where the Mach number is equal to unity, and the equation D(r, M, C) =0 determines its location. 

The extra regularity conditions at the sonic point Afi(r,Ai,C) — are satisfied only for one 
particular value of the angular momentum at the horizon £j n , which is the eigenvalue of the 
problem that should be found. For a given a the location of the sonic point depends on the 
mass accretion rate. For low mass accretion rates one expects the transonic transition to occur 
close to the ISCO. Figure [5] shows that this is indeed the case for accretion rates smaller than 
about 0.4AfEdd, independent of a, where we use the authors' definition of AfEdd = 16LEdd/c 2 - 
At M = 0.4MEdd a qualitative change occurs, resembling a "phase transition" from the Shakura- 
Sunyaev behavior to a very different slim-disk behavior. For higher accretion rates the location 
of the sonic point significantly departs from the ISCO. For low values of a, the sonic point moves 
closer to the horizon, down to ~ 4M for a = 0.001. For a > 0.2 the sonic point moves outward 
with increasing accretion rate, reaching values as high as 8M for a = 0.5 and M = 100 AfEdd' 
This effect was first noticed for small accretion rates by Muchotrzeb |212j and later investigated 
for a wide range of accretion rates by Abramowicz |10j , who explained it in terms of the disk-Bondi 
dichotomy. 

The topology of the sonic point is important, because physically acceptable solutions must be 
of the saddle or nodal type; the spiral type is forbidden. The topology may be classified by the 
eigenvalues Ai, A2, A3 of the Jacobi matrix, 



J 



dV dV &D 

dr, dM dC 

aVi 0N\ dNi 

dr. dM dC 

aVa M2 cWk 

dr dM dC 



(95) 



Because det(^7) = 0, only two eigenvalues Ai,A2 are non-zero, and the quadratic characteristic 
equation that determines them takes the form, 

2A 2 -2Atr(,7)- [tr(.7 2 )-tr 2 (J)] =0. (96) 
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Figure 5: Location of the sonic point as a function of the accretion rate for different values of a, for 
a non-rotating black hole, a — 0, taking MEdd = 16^EddC 2 . The solid curves are for saddle type 
solutions while the dotted curves present nodal type regimes. Image reproduced by permission 
from [9], copyright by ESO. 



The nodal-type solution is given by A1A2 > and the saddle type by A1A2 < 0, as marked in 
Figure [5] with the dotted and the solid lines, respectively. For the lowest values of a only the 
saddle-type solutions exist. For moderate values of a (0.1 < a < 0.4) the topological type of the 
sonic point changes at least once with increasing accretion rate. For the highest a solutions, only 
nodal- type critical points exist. 



5.3 Solutions: Shakura Sunyaev &z Novikov— Thorne 



Shakura and Sunyaev [279 noticed that a few physically reasonable extra assumptions reduce the 
system of thin disk equations (88 1 -(93) to a set of algebraic equations. Indeed, the continuity 



and vertical equilibrium equations, (88) and (92), are already algebraic. The radial momentum 



equation ( 90 ) becomes a trivial identity = with the extra assumptions that the radial pressure 
and velocity gradients vanish, and the rotation is Keplerian, f2 = Q^. The algebraic angular 
momentum equation (91) only requires that we specify Ci n . The Shakura-Sunyaev model makes 
the assumption that £j„ = £fc(ISCO). This is equivalent to assuming that the torque vanishes 
at the ISCO. This is a point of great interest that has been challenged repeatedly [1641 11041 125] . 



Direct testing of this hypothesis by numerical simulations is discussed in Section |11.4 



The right-hand side of the energy equation ( 93 1 represents advective cooling. This is assumed 



to vanish in the Shakura-Sunyaev model, though we will see that it plays a critical role in slim disks 
(Section [6]) and ADAFs (Section W\. Because the Shakura-Sunyaev model assumes the rotation 
is Keplerian, f2 = Qt , meaning i7 is a known function of r, the first term on the left-hand side 



of Eq. (93), which represents viscous heating, is algebraic. The second term, which represents the 



radiative cooling (in the diffusive approximation) is also algebraic in the Shakura-Sunyaev model. 

In addition to being algebraic, these thin-disk equations are also linear in three distinct radial 
ranges: outer, middle, and inner. Therefore, as Shakura and Sunyaev realized, the model may be 
given in terms of explicit algebraic (polynomial) formulae. This was an achievement of remarkable 
consequences - still today the understanding of accretion disk theory is in its major part based 
on the Shakura-Sunyaev analytic model. The Shakura-Sunyaev paper [279] is one of the most 
cited in astrophysics today (see Figure pj), illustrating how fundamentally important accretion disk 
theory is in the field. 

The general relativistic version of the Shakura-Sunyaev disk model was worked out by Novikov 
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Figure 6: The number of citations to the Shakura & Sunyaev paper |279j is still growing exponen- 
tially, implying that the field of black hole accretion disk theory still has not reached saturation. 
Image reproduced from the SAO/NASA Astrophysics Data System: http://adsabs. harvard. 
|edu/abs/ 1973AM 24. .337S[ 



and Thorne [229) . with important extensions and corrections provided in subsequent papers [2371 
265, 1241) . Here we reproduce the solution, although with a more general scaling: m = M/M© and 
m = Mc 2 /L Edd . 

Outer region: P = P gas , k — Kff (free-free opacity) 



F 
£ 
H 

Po 
T 

0/(1-0) 

T ff/ T es 



[7 x 10 26 erg cm" 2 s _1 ](m _1 m) r~ 3 S" 1 C" 1/2 Q, 

[4 x 10 5 g cm- 2 ]( a - 4 / 5 m 2 / 10 m^ 10 )n 3/4 ^ 1 /io s -4/5 C i/2 I3 -i7/20 f -i/2o g 7/iO ; 
[4 x 10 2 cm]^- 1 / 10 ™ 18 / 20 ™ 3 / 2 ^ 7 ^^ 

[4 x 10 2 g cm- 3 ]^" 7 / 10 ™- 7 / 10 ™ 11 / 20 )^ 15/8 ^- 17/20 B 3/1 °p- 11/40 f 17/40 Q 11/20 
[2 x 10 s K](a~ 1/5 m~ 1/5 7h 3/1 °) r - 3 / 4 A- 1/10 B~ 1/a 'D- 3/20 £ 1/2a Q 3/10 , 
[3}(a-^ 10 m-^ 10 m- 7 / 20 )rl /8 ^-n^/io^o^ ii/40 g -7/2 0j 
[2 x 10- 3 ](m- 1 / 2 )r 3/4 ^- 1 /2 6 2/5 23 i/4 f i/4 g -i/2 ; 



(97) 



where = rc 2 /GM. 

Middle region: P = P S3 , S , K = n es (electron-scattering opacity) 



F 
£ 
H 

Po 
T 

T ff/ T es 



[7 x 10 26 erg cm" 2 s^Km^m) r~ 3 S^C^S, 

[9 x 10 4 g cm- 2 ](a- 4 / 5 m 1 / 5 m 3 / 5 )r- 3/5 6- 4 / 5 C 1 / 2 2?- 4 / 5 Q 3 / 5 , 

[1 x 10 3 cm](a- 1 / 10 m 9 / 10 m 1 / 5 )r 21/20 ^- 6 / 5 C 1 / 2 2?- 3 / 5 f- 1 / 2 Q 1 / 5 , 
[4 x 10 1 g cm- 3 ](a- 7 / 10 m- 7 / 10 ™ 2 / 5 )n 33/20 ^- 1 S 3 / 5 2?- 1 / 5 £ 1 /2 S 2/5 ! 

[7 x 10 8 K](a _1/5 m" 1/5 m 2/5 ) rZ 9/1 ° B _2 / 5 D- 1 / 5 Q 2 / 5 ) 

[7 x 10- 3 ](«- 1 / 1 V- 1 / 10 m- 4 / 5 )r 21/20 ^- 1 6 9 / 5 2? 2 / 5 £ 1 / 2 Q- 4 / 5 , 

[2 x W^m- 1 ) r 3 / 2 A- 1 B 2 V 1 / 2 £ 1 / 2 Q-\ (98) 
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Inner region: P = P ra( j, K = K es 



F 
E 
H 

Po 
T 

(3/(1- (3) 
(r//r e ,) 1/2 



[7 x 10 26 erg cm~ 2 s^Km^m) r~ 3 B^C^^Q, 
[5 g cm" 2 ]^ 1 ™- 1 ) r 3 J 2 A- 2 B 3 C^ 2 £Q-\ 
[1 x 10 5 cm}(m)A 2 B- 3 C 1/2 V- 1 £- 1 Q, 



= [2 x 10" 



g cm ](a ' /// ' in 



= [5 X 10 7 K](o- 1 /* m -V4) r -3/8 i/r l/2 jB l/2 f l/4 ) 

= [4 x lO-^la- 1 /^- 1 / 4 ™- 2 )^ 178 ^- 5 / 2 ^/ 2 ^ 5 /^- 2 , 

= [1 x 10- 4 ](a- 17 / 16 m^ 1 / 16 m- 2 )rf /32 ^ 25 / 8 6 41 / 8 C 1 / 2 I? 1 / 2 £ 25 / 16 Q- 2 . (99) 



The radial functions .4, Q that appear in Eqs. (97), (98), (99), are defined in terms of y 
(r/M) 1 / 2 and a* = a/M as [237]: 



2 -4 

a *2/ 



C = 1 - 3y~ 2 + 2a*y -3 , 
£ = l + 4aly- 4 - 4a^- 6 



B = 

D = 

Qo = 



-2 



2y" 



2 -4 

a*2/ , 



1 + a*y~ 



2 /(l-3y- 2 +2a,?/- 3 ) 1 /2' 
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V2/0 


-2/3/. 



IS 



Here ?/o = (^ms/-^) 1 , and j/i, 2/2, an d j/ 3 are the three roots of y 3 — 3y + 2a* = 0; that 

yi = 2cos[(cos~ 1 a* — 7r)/3], 
y 2 = 2cos[(cos _1 a* + 7r)/3], 
2/3 = -2cos[(cos^ 1 a*)/3]. 

For the numerical solutions reproduced in Figure [7| the opacities were assumed to be n es = 
0.34 cm 2 g _1 and Kft = 6.4 x 10 22 /9 cgs Tj^ 7 ^ 2 cm 2 g -1 , where p cgs is the density in g cm -3 and 7k 
is the temperature in Kelvin. 

The Shakura-Sunyaev and Novikov-Thorne solutions are only local solutions; this is because 
they do not take into account the full eigenvalue problem described in Section [572} Instead, they 
make an assumption that the viscous torque goes to zero at the ISCO, which makes the model 
singular there. For very low accretion rates, this singularity of the model does not influence the 
electromagnetic spectrum [298] . nor several other important astrophysical predictions of the model. 
However, in those astrophysical applications in which the inner boundary condition is important 
(e.g., global modes of disk oscillations), the Novikov-Thorne model is inadequate. Figure [7] illus- 
trates a few ways in which the model fails to capture the true physics near the ISCO. 
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Figure 7: The innermost part of the disk. In the Shakura-Sunyaev and Novikov-Thorne models, 
the locations of the maximum pressure (a.k.a. the center) r con t or and the cusp r cusp , as well as the 
sonic radius r soun( j, are assumed to coincide with the ISCO. Furthermore, the angular momentum 
is assumed to be strictly Keplerian outside the ISCO and constant inside it. In real flows, r centor ^ 
^cusp 7^ ''sound 7^ ISCO, and angular momentum is super-Keplerian between r cusp and r ccn tcr- 
Image reproduced by permission from [5], copyright by ESO. 
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6 Slim Disks 



The Shakura-Sunyaev and Novikov-Thorne models of thin disks assume that accretion is radia- 
tively efficient. This assumption means that all the heat generated by viscosity at a given radius 
is immediately radiated away. In other words, the viscous heating is balanced by the radiative 
cooling locally and no other cooling mechanism is needed. This assumption can be satisfied as 
long as the accretion rate is small. At some luminosity (L f=a 0.3 ^Edd), however, the radial velocity 
is large, and the disk is thick enough, to trigger another mechanism of cooling: advection. It 
results from the fact that the viscosity-generated heat does not have sufficient time to transform 
into photons and leave the disk before being carried inwards by the motion of the gas. The higher 
the luminosity, the more significant advective cooling is. At the highest luminosities, it becomes 
comparable to the radiative cooling (see Figure [8]) , and the standard, thin disk approach can no 
longer be applied. 




Figure 8: The advection factor (ratio of advective to radiative cooling) profiles for m = 0.01, 1.0 
and 10.0 (here, m = M c 2 / '16 L^d) ■ Profiles for a = 0.01 and 0.1 are presented with solid black 
and dashed red lines, respectively. The fraction / a dv/(l + /adv) of heat generated by viscosity 
is carried along with the flow. In regions with / adv < the advected heat is released. Image 
reproduced by permission from [269) . 

The problem of accretion with an additional cooling mechanism has to be treated in a different 
way than radiatively efficient flows. Without the assumptions of radiative efficiency and Keplerian 
angular momentum, it is no longer possible to find an analytic solution to the system of equations 
presented in Section [5.1| Instead, one has to solve a two-dimensional system of ordinary differential 
equations ( |94| with a critical point - the radius at which the gas velocity exceeds the local speed 
of sound (the sonic radius). This was first done in the pseudo-Newtonian limit by Abramowicz [5], 
who forged the term "slim disks" . It has since been done using a fully relativistic treatment by 
Bcloborodov |40J. Recently, Sadowski 268 constructed slim disk solutions for a wide range of 
parameters applicable to X-ray binaries. 

These slim disks are in some sense more physical than thin disks and offer a more general 
set of solutions, while in the limit of low accretion rates they converge to the standard thin disk 
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solutions. Slim disks are more physical in that they extend down to the black hole horizon, as 
opposed to thin disks that formally terminate at the ISCO. Slim disks are more general in that 
they may rotate with an angular momentum profile significantly different than the Keplerian 
one - the higher the accretion rate, the more significant the departure (see Figure [9]). The disk 
thickness also increases with the accretion rate. For rates close the Eddington limit, the maximal 
H/R ratio reaches 0.3. Finally, the flux emerging from the slim disk surface is modified by the 
advection. At high luminosities, a large fraction of the viscosity-generated heat is advected inward 
and released closer to the black hole or not released at all. As a result, the slope of the radial flux 
profile changes, and radiation is even emitted from within the ISCO (see Figure 10 1. Due to the 
increasing rate of advection, the efficiency of transforming gravitational energy into radiative flux 
decreases with increasing accretion rate. Despite highly super-Eddington accretion rates, the disk 
luminosity remains only moderately super-Eddington (see Figure 11). The Eddington luminosity 
may be exceeded because the geometry of the flow is not spherical and the classical definition of 
this quantity does not apply - most of the accretion takes place in the equatorial plane while the 
radiation escapes vertically. Thus, the radiation is not capable of stopping the inflow, though it 
may cause outflows from the surface. 
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Figure 9: Profiles of the disk angular momentum (u^) for a = 0.01 (left) and a = 0.1 (right panel) 
for different accretion rates (as a reminder, m = M c 2 / 16 L^d) , showing the departures from the 
Keplerian profile. These plots are for a non-rotating black hole. Image reproduced by permission 
from [270] . copyright by ESO. 
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Figure 10: Flux profiles for different mass accretion rates in the case of a non-rotating black hole 
and two values of a: 0.01 (black solid), 0.1 (red dashed lines). For each value of a there are five 
lines corresponding to the following mass accretion rates: rh = 0.01, 0.1, 1.0, 2.0 and 10.0 (as a 
reminder, m = Afc 2 /16 ^Edd)- The black hole mass is 10 M Q . Image reproduced by permission 
from [UD], copyright by ESO. 
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Figure 11: Top panel: Luminosity vs accretion rate for three values of black hole spin (a* = a/M = 
0.0, 0.9, 0.999) and two values of a = 0.01 (black) and 0.1 (red line). Bottom panel: efficiency of 
accretion r/ = (L / '^Edd) / (M / 'MEdd) (here M^dd — 16£Edd/c 2 )- Irnage reproduced by permission 
from [259]. 
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7 Advection-Dominated Accretion Flows (ADAFs) 



The ADAF, or advection-dominated accretion flow, solution also involves advective cooling. In 
fact, it carries it to an extreme - nearly all of the viscously dissipated energy is advected into the 
black hole rather than radiated. Unlike the slim disk solution, which is usually invoked at high 
luminosities, the ADAF applies when the luminosity (and generally the mass accretion rate) are 
low. 

Because of their low efficiency, ADAFs are much less luminous than the Shakura-Sunyaev thin 
disks. The solutions tend to be hot (close to the virial temperature), optically thin, and quasi- 
spherical (see Figure 12 1. Their spectra are non-thermal, appearing as a power-law, often with a 
strong Compton component. This makes them a good candidate for the Hard state observed in 



X-ray binaries (discussed in Section 12.3 ) 
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Figure 12: Profiles of temperature, optical depth, ratio of scale height to radius, and advection 
factor (the ratio of advective cooling to turbulent heating) of a hot, one-T ADAF (solid lines). 
The parameters are M = 10 M Q , M = 1(T 5 L Edd /c 2 , a = 0.3, and /3 = P gas /(P glls + P mag ) = 0.9. 
The outer boundary conditions are R out — 10 3 Rs, T — 10 9 K, and v/c s — 0.5. Two-T solutions 
with the same parameters and S = 0.5 (dashed lines) and 0.01 (dot-dashed lines) are also shown 
for comparison, where S is the fraction of the turbulent viscous energy that directly heats the 
electrons. Image reproduced by permission from [321] . copyright by AAS. 

ADAFs were formally introduced in the Newtonian limit through a series of papers by Narayan 
and Yi [223, 224, 225 , followed closely by Abramowicz [UE] and others |106j . although the existence 
of this solution had been hinted at much earlier [134[ I257j . In the same spirit as we gave the 
equations for the Novikov-Thorne solution in Section [5.3| for thin disks, we report the self-similar 
ADAF solution found by Narayan and Yi [224] . Again we present the solution with the following 
scaling: r* = rc 2 /GM, m = M/M Q and m = Mc 2 /L E dd- 
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where i> is the radial infall velocity and q + is the viscous dissipation of energy per unit volume. 
The constants c\, C2 , and C3 are given by 
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and the parameter / a d v represents the fraction of viscously dissipated energy which is advected. 
The remaining amount, 1 — / a dv, is radiated locally. 

The rapid advection in ADAFs generally has two effects: 1) dissipated orbital energy can not 
be radiated locally before it is carried inward and 2) the rotation profile is generally no longer 
Kcplerian, although Abramowicz [5] found solutions where the dominant cooling mechanism was 
advection, even when the angular momentum profile was Kcplerian. Fully relativistic solutions of 
ADAFs have also been found numerically [T31 HI] . Further discussion of ADAFs is given in the 
review article by Narayan and McClintock |219j . 
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8 Stability 



Having reviewed some of the main analytic models of accretion disks, it is important now to discuss 
the issue of stability. Since all analytic models presume steady-state solutions, such models are 
only useful if the resulting solutions are stable. One reason to suspect accretion disks may not 
be stable is that the systematic differential rotation that is a signature feature of accretion is a 
potential source of energy, and therefore, of instability. Another is that some level of instability 
may be essential in accretion disks as it can provide a pathway to the kind of sustained turbulence 



anticipated by Shakura and Sunyaev (see Section 3.2.11 



8.1 Hydrodynamic stability 

Within ideal hydrodynamics, local linear stability of an axisymmetric rotating flow is guaranteed 
if the H0iland criterion is satisfied 



1 di 2 1 

WdR~C^ p - VS > °' (1 ° la) 



d P (d£ 2 dS d£ 2 dS\ 



< 0, (101b) 



where C v is the specific heat at constant pressure and R is the cylindrical radius (see |275j for 
the criterion for relativistic stars). This criterion can be easily understood in two limits: For 
non-rotating equilibria (e.g., a non-rotating star), the criterion reduces to the Schwarzschild cri- 
terion {dS/dr > 0) that the entropy must not increase toward the interior (for stability against 
convection). Provided this is true, local fluid elements will simply oscillate under stable buoyancy 
forces. To see the effects of rotation, we can consider an equilibrium that has constant entropy 
everywhere. Then the H0iland criterion reduces to the Rayleigh criterion (d£ 2 /dR > 0): the spe- 
cific angular momentum must not decrease outward. Physically, if one perturbs a fluid element 
radially outward, it conserves its own specific angular momentum. If the ambient specific angular 
momentum decreases outward, then the fluid element will be rotating too fast to stay in its new 
position, and centrifugal forces will push it further outward. Stability would be a fluid element 
that oscillates at the local epicyclic frequency. 

As it turns out, the H0iland criterion is a huge disappointment for understanding why turbu- 
lence might exist in accretion disks. This is because it indicates that accretion disks with rotation 
profiles that do not differ too much from Keplerian should be strongly stable! 



8.1.1 Papaloizou Pringle Instability (PPI) 

The H0iland criterion is only a local stability criterion. Flows can be locally stable, yet have global 
instabilities. An example of this occurs in the Polish doughnut solution (Section [4]) . Papaloizou 
and Pringle [238j showed that this solution is marginally stable with respect to local axisymmetric 
perturbations yet unstable to low-order nonaxisymmetric modes. As with all global instabilities, 
the existence of the Papaloizou-Pringle instability (PPI) is sensitive to the assumed boundary 
conditions [44]. In cases where the disk overflows its potential barrier (Roche lobe) and accretes 
through pressure-gradient forces across the cusp, the PPI is generally suppressed |117j . 



8.1.2 Runaway instability 

Another instability associated with the Polish doughnut is the runaway instability [5] . If matter is 
overflowing its Roche lobe and accreting onto the black hole, then one of two evolutionary tracks 
are possible: (i) As the disk loses material it contracts inside its Roche lobe, slowing the mass 
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transfer and resulting in a stable situation, or (ii) as the black hole mass grows, the cusp moves 
deeper inside the disk, causing the mass transfer to speed up, leading to the runaway instability. 
Recent numerical simulations show that, while this instability grows very fast, on timescales of a 
few orbital periods, over a wide range of disk-to-black hole mass ratios when i = const., i.e., a 
constant specific angular momentum profile 98 , it is strongly suppressed whenever the specific 
angular momentum of the disk increases with the radial distance as a power law, I oc r p [63] . Even 
values of p much smaller than the Keplcrian limit (p = 1/2) suffice to suppress this particular 
instability. [This is equivalent to angular velocity profiles, SI oc r~ q , with q > 3/2.] 

8.2 Magneto-rotational instability (MRI) 

Although it had long been suspected that some sort of MHD instability might provide the necessary 
turbulent stresses to make accretion work, the nature of this instability remained a mystery until 
the rediscovery of the magneto-rotational instability by Balbus and Hawley [26llll8|, l2"7]. Originally 
discovered by Velikhov |309) . and generalized by Chandrasekhar [SB], in the context of vertically 
magnetized Couette flow between differentially rotating cylinders, the application of this instability 
to accretion disks was originally missed. 

The instability itself can be understood through a simple mechanical model. Consider two 
particles of gas connected by a magnetic field line. Arrange the particles such that they are initially 
located at the same cylindrical distance from the black hole but with some vertical separation. 
Give one of the particles (say the upper one) a small amount of extra angular momentum, while 
simultaneously taking away a small amount of angular momentum from the lower one. The upper 
particle now has too much angular momentum to stay where it is and moves outward to a new 
radius. The lower particle experiences the opposite behavior and moves to a smaller radius. In 
the usual case where the angular velocity of the flow drops off with radius, the upper particle will 
now be orbiting slower than the lower one. Since these two particles are connected by a magnetic 
field line, the differing orbital speeds mean the field line will get stretched. The additional tension 
coming from the stretching of the field line provides a torque, which transfers angular momentum 
from the lower particle to the upper one. This just reinforces the initial perturbation, so the 
separation grows and angular momentum transfer is enhanced. This is the fundamental nature of 
the instability. 

In more concrete terms, consider a disk threaded with a vertical magnetic field B z and having 
an Alfven speed v\ = B 2 /(4irp). The dispersion relation for perturbations of a fluid quantity 
SX ~ exp[i(kz - u)t)] is [27] 

w 4 - (2kv A + uj 2 )uj 2 + kv A (kv A + rdtf/dr) = . (102) 

This equation has an unstable solution (lu 2 < 0), if and only if, kv A + rdfl 2 /dr < 0. Since 

<9f2 

t < °- ^ 
in accretion disks, the instability criteria can generally be met for weakly magnetized disks. More 
specifically, the MRI exists for intermediate magnetic field strengths. In terms of the natural 
length scale of the instability (~ v A /Q), the field strength is constrained at the upper limit by 
the requirement that the unstable MRI wavelength fit inside the vertical thickness of the disk 
(v A /Q < H). This corresponds to field energy densities that are less than the thermal pressure, 
i.e., b 2 < -Pgas- At the lower end, dissipative processes set a floor on the relevant length scales, and 
hence, field strengths. 

If the conditions for the instability are met, the fastest-growing mode, which dominates the early 
evolution, has the form of a "channel flow" involving alternating layers of inward- and outward- 
moving fluid. The amplitude of this solution grows exponentially until it becomes unstable to three- 
dimensional "parasitic modes" that feed off the gradients of velocity and magnetic field provided by 
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the channel flow. The flow rapidly reaches a state of magnetohydrodynamic turbulence |1181 1119] . 
This instability can be self-sustaining through a nonlinear dynamo process [55] - nonlinear because 
the motion that sustains or amplifies the magnetic field is driven by the field itself through the 
MRI. A more complete description of the linear and non-linear evolution of the MRI is provided 
in the review article by Balbus and Hawley [27]. A general relativistic linear analysis is presented 
in [2UJ. 



8.3 Thermal and viscous instability 

It was realized by Shakura and Sunyaev themselves [280] . as well as other authors |176l 1287] . 



that the Shakura-Sunyaev solution (Section 5.3) should be thermally and viscously unstable for 



disks in which radiation pressure dominates (when the opacity is governed by electron scattering) . 
The most general and elegant arguments are presented in the classic paper by Piran |246j . This 
discovery started a long debate, which continues unresolved to this day. A recent update is provided 
in [50]. 

To understand the thermal instability better, we consider a disk cooling through radiative 
diffusion. The local emergent flux at radius r is given by 

F~ = ^, (104) 

T 

where a is the radiation density constant, T is a measure of the disk interior temperature, and r 
is half the total vertical optical depth. If the opacity is dominated by electron scattering, then 
r ~ k^S/2, where £ is the surface density of the disk, which is constant on the time scales of 
the instabilities (very much less than the radial flow time scale). The Thomson opacity kt is 
also constant, being independent of temperature, provided there is already sufficient ionization. 
Therefore, the optical depth is independent of temperature and the cooling rate per unit area is 
F~ oc T 4 . 

The dissipation rate per unit area is 

F + ~rHT r ^. (105) 
Vertical hydrostatic equilibrium implies that the disk half thickness H ~ 2P/(f2| f £), so that 



Eq. (1051 becomes 



In the radiation pressure dominated inner region, P ~ aT /3, so that Eq. (106) plus the Shakura- 
Sunyaev assumption T r $ = —aP imply that F + oc T 8 l Hence a perturbative increase in tempera- 
ture increases both the local cooling and heating rates, but the heating rate increases much faster, 
leading to a thermal runaway. 

Note, though, that this argument only applies when the viscous stress is proportional to the total 
pressure Prot (p. being the proportionality constant). For some time it seemed that a plausible way 
to avoid this instability was to argue that the stress is proportional instead to the gas pressure P gas . 
Recent numerical simulations, though, of the magneto-rotational instability in radiation-pressure 
dominated disks have shown that the stress is, in fact, proportional to the total pressure [128] . 
Interestingly, these simulations exhibit no sign of the predicted thermal instability. 

Most observations also argue against the existence of this instability. In the case of accretion 
onto black holes, the instability is supposed to set in for luminosities in excess of L > 0.01 pEdd- 
However, during outbursts, many stellar-mass black hole sources cross this limit both during their 
rise to peak luminosity and on their decline to quiescence, showing no dramatic symptoms (although 
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they do undergo state changes, as described in Section 12.3). On the contrary, observations suggest 
that disks in black hole X-ray binaries are stable up to at least L ss 0.5LEdd [HS]- Certainly there 
is no evidence for the sensational behavior anticipated by some models |172l 1300] . 
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9 Oscillations 



Even when analytic disk solutions are stable against finite perturbations, it is often the case that 
these perturbations will, nevertheless, excite oscillatory behavior. Oscillations are a common dy- 
namical response in many fluid (and solid) bodies. Here we briefly explore the nature of oscillations 
in accretion disks. This topic is particularly relevant to understanding the physical mechanisms 
that may be behind quasi-periodic oscillations (or QPOs, which are discussed in Section 12.4). 

There are a number of local restoring forces available in accretion disks to drive oscillations. 
Local pressure gradients can drive oscillations via sound waves. Buoyancy forces can act through 
gravity waves. The Coriolis force can operate through inertial waves. Surface waves can also exist, 
with the restoring force given by the local effective gravity. 

Of particular interest are families of low order modes that may exist in various accretion ge- 
ometries. Such modes will tend to have the largest amplitudes and produce more easily observed 
changes than their higher-order counterparts. Here we briefly review a couple relatively simple 
examples for the purpose of illustration. More details can be found in the references given. 



9.1 Dynamical oscillations of thick disks 

A complete analysis of the spectrum of modes in thick disks has not yet been done. Some progress 
has been made by considering the limiting case of a slender torus, where slender here means 
that the thickness of the torus is small compared to its radial separation from the central mass 
(i.e. the torus has a small cross-sectional area). In this limit, the complete set of modes have 
been determined for the case of constant specific angular momentum in a Newtonian gravitational 
potential . A more general analysis of slender torus modes is given in [15J . 

Any finite, hydrodynamic flow orbiting a black hole, such as the Polish doughnuts described in 
Section|4j is susceptible to axisymmetric, incompressible modes corresponding to global oscillations 
at the radial (a — uj r ) and vertical (cr = luq) epicyclic frequencies. Other accessible modes are found 
by solving the relativistic Papaloizou-Pringle equation [3] 
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m is the azimuthal wave number, n is the polytropic index (assuming an equation of state of the 
form P = Kp 1+1 / n ), c s o is the sound speed at the pressure maximum r , and g is the determinant of 
the metric. The necessary boundary condition comes from requiring that the Lagrangian pressure 
perturbation vanish at the unperturbed surface (/ = 0): 



AP = (6P + £ a V a P) = 0. 



(109) 



where is the Lagrangian displacement vector. Eqs. (1071 and (109) describe a global eigen- 



value problem for the modes of Polish doughnuts, with three characteristic frequencies: the radial 
epicyclic frequency cj r , the vertical epicyclic frequency loq, and the characteristic frequency of 
inertial modes k. 
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In the slender-torus limit, one can write down analytic expressions for a few of the lowest 
order modes [25] besides just the radial and vertical epicyclic ones. In terms of local coordinates 
measured from the equilibrium point, 

an eigenfunction of the form W = Axy, for some constant A, yields two modes with eigenfrequen- 
cies [45] 



\ [u 2 , + u 2 ± [K 2 + u 2 ) 2 + 4,&f\ 1/2 } , (111) 

where ao = ctq/CIq. The positive square root gives a surface gravity mode that has the appearance 
of a cross(x )-mode. The negative square root gives a purely incompressible inertial (c-) mode, 
whose poloidal velocity field represents a circulation around the pressure maximum. 



Radial Epicyclic (- + 01) Vertical Epicyclic ( + -01) 




x Mode (--02) Inertial Mode ( — 02) 




Figure 13: Poloidal velocity fields (Su x , 8u y ) of the lowest order, non-trivial thick disk modes. 
Image reproduced by permission from [45 j . 

An eigenfunction of the form W = A + Bx 2 + Cy 2 , with arbitrary constants A, B, C, also has 
two modes with eigenfrequencies |45j 

a 2 = i{( 2 » + l)( w 3+ w 2)-(n+l)«g (112) 
±[((2n + 1){lo 2 - u 2 r f + (n + l)4) 2 + 4(u r 2 - $)u$\V*} . 

The upper sign results in an acoustic mode with the velocity field of a breathing mode. This mode 
is comparable to the acoustic mode in the incompressible Newtonian limit for I = const., while in 
the Keplerian limit, the mode frequency becomes that of a vertical acoustic wave. The lower sign 
corresponds to a gravity wave, with a velocity field reminiscent of a plus(+)-mode. The poloidal 
velocity fields of all these lowest order modes are illustrated in Figure [l3j 

9.2 Diskoseismology: oscillations of thin disks 

To analyze oscillations of thin disks, one can express the Eulerian perturbations of all physical quan- 
tities through a single function SW cx SP/p. Since the accretion disk is considered to be stationary 
and axisymmetric, the angular and time dependence are factored out as SW — W(r, z)e l ^ m ^ crt \ 
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where the eigenfrequency a{r, z) = oj — mfl. Here it is useful to assume that the variation of the 
modes in the radial direction is much stronger than in the vertical direction, z = r cos 9. The re- 
sulting two (separated) differential equations for the functional amplitude W(r, z) — W r (r)W y (r, y) 
are given by [221 EH E2] 



d 2 W r 

dr 2 (co 2 



dW r , (« t ?9r* [ J_ uT] ( 1 _ *W =0 , (113 ) 



dr c s r \ U) 



(1 -V 2 )^r~ ^y^f + 2g [* 2 y 2 + - y 2 )] W y = 0. (114) 

The radial eigenf unction, W r , varies rapidly with r, while the vertical eigenfunction, W y , varies 
slowly. Here y — [z /H)\y g / (jg — l)] 1 ^ 2 is the re-scaled vertical coordinate, 7 9 is the adiabatic 



index, (I>(r) = Lu r /ujg is the ratio of the epicyclic frequencies from Section 2.3 $ is the eigenvalue 
of the (WKB) separation function, and g = (P/P c )/(p/p c ), where P c and p c are the midplane 
values. The radial boundary conditions depend on the type of mode and its capture zone (see 
below). Oscillations in thin accretion disks, then, are described in terms of ^(r, a), along with the 
angular, vertical, and radial mode numbers (number of nodes in the corresponding eigenfunction) 
m, j, and n, respectively. Modes oscillate in the radial range where 

(" 2 -^)(l- J) >0 (115) 

outside the inner radius, r > r\. 

p-modes are inertial acoustic modes defined by $ < Cj 2 and are trapped where oj 2 > oj 2 , 
which occurs in two zones. The inner p-modes are trapped between the inner disk edge and the 
inner "Lindblad" radius, i.e., n < r < r_, where gas is accreted rapidly. The outer p-modes occur 
between the outer Lindblad radius and the outer edge of the disk, i.e., r + < r < r . The Lindblad 
radii, r_ and r + , occur where uj = u r . The outer p-modes are thought to be more consequential as 
they produce stronger luminosity modulations [233) . In the corotating frame these modes appear 
at frequencies slightly higher than the radial epicyclic frequency. Pressure is the main restoring 
force of p-modes. 

g- modes are inertial gravity modes defined by W > uj 2 . They are trapped where u> 2 < ui 2 
in the zone r_ < r < r + given by the radial dependence of ui ri i.e., g- modes are gravitationally 
captured in the cavity of the radial epicyclic frequency and are thus the most robust among the 
thin-disk modes. Since this is the region where the temperature of the disk peaks, g-modes are also 
expected to be most important observationally [244]. In the corotating frame these modes appear 
at low frequencies. Gravity is their main restoring force. 

c-modes are corrugation modes defined by W = ui 2 . They are non-radial (m = 1) and vertically 
incompressible modes that appear near the inner disk edge and precess slowly around the rotational 
axis. These modes are controlled by the radial dependence of the vertical epicyclic frequency. In 
the corotating frame they appear at the highest frequencies. 

All modes have frequencies cx 1/M. Upon the introduction of a small viscosity (v <x a, a 1), 
most of the modes grow on a dynamical timescale idym such that the disk should become unstable. 
However, evidence for these modes has so far mostly been lacking in MRI turbulent simulations 



(see Section 11.6). This leaves their relevance in some doubt. 
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10 Relativistic Jets 



Although the main focus of this review is on black hole accretion disk theory, we note that there 
has long been a strong observational connection between accreting black holes and relativistic jets 
across all scales of black hole mass. For supermassive black holes this includes quasars and active 
galactic nuclei; for stellar-mass black holes this includes microquasars. However, the theoretical 
understanding of disks and jets has largely proceeded separately and the physical link between the 
two still remains uncertain. Therefore, we present only a few brief comments on the subject in 
this review. More complete discussions of the theory of relativistic jets may be found in |193| . A 
review of their observational connection to black holes is given in [205 . 



In Section 2.2 we described the "Penrose process," whereby rotational energy may be extracted 
from a black hole and carried to an observer at infinity. To briefly recap, Penrose [242. 243 imagined 
a freely falling particle with energy E°° disintegrating into two particles with energies 25f° < and 
ET > 0. Then, the particle with negative energy Ef falls into the black hole, and the other one 
escapes to infinity. Clearly, > E°° , so that there is a net gain of energy. 

It was first suggested by Wheeler at a 1970 Vatican conference and soon after by others [1831 
I96j that such a Penrose process may explain the energetics of superluminal jets commonly seen 
emerging from quasars and other black hole sources. However, a number of authors [31, 314, 161 
showed that for to be greater than E°°, the disintegration process must convert most of the 
rest mass energy of the infalling particle to kinetic energy, in the sense that, in the center-of-mass 
frame, the E^ particle must have velocity v > c/2. The argument of Wald [314J is powerful, short 
and elegant, so we give it here in extenso. 

Let p v = mu v be the four momentum of a particle with the mass m. We assume that in the 
ZAMO frame (Section |2~2| the particle has a four velocity of the form, u v = T(N V + VA U ), where 
A" is a timelike- unit vector (for simplicity we assume A" r\ v = 0), V is the particle 3- velocity in 
the ZAMO frame, and T~ 2 = 1 — V 2 . If the disintegration fragments move in the directions ±A W 
(which one may prove is energetically most favorable), then the four velocities of these fragments 
in the center-of-mass frame are, 

u v ± = 1 {u v ±vt v ), (116) 
t v = T(VN V +A U ). (117) 



The form of (117) follows from the Lorentz transformation {N u ,k u } — > {u",t"}. Multiplying 
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Since l/g u < 1 and realistic particles have E/m > l/v3> the condition E°° < necessarily 
requires v > 1/2. Such highly relativistic disintegration events are not generally seen in nature. To 



make matters worse, from the upper limit of (118), Wald deduced that the presence of the black 



hole limits the energy increase to a maximal factor of 1 + y/2. Thus, he concluded |314j : "The 
Penrose mechanism cannot serve as a useful energy source for astrophysical processes. In no case 
can one obtain energies which are greater by a significant factor than those which already could be 
obtained by a similar breakup process without the presence of the black hole." 

Replacing particle disintegration with particle collision does not help, even though the center-of- 
mass energy of such a collision happening arbitrarily close to the horizon of the maximally rotating 
Kerr black hole may be arbitrarily large |247[ [28] . This is because the Wald limit of 1 + v2 still 
holds |39j . It would seem that even under idealized conditions, the maximal energy of a particle 
escaping via the Penrose process is only a modest factor above the total initial energy [33] . 

Therefore, we consider a general matter distribution, described by an unspecified stress-energy 
tensor Ttf. In this case, the energy flux in the ZAMO frame is E % = —T l k N k , and the energy 
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absorbed by the black hole is 

E= - I T k n K dNi > 0, (119) 



where J dNi is the surface integral over the horizon. The inequality sign follows from the fact that 
the locally measured energy must be positive. The above integral may by transformed into 

< E = - J e* T\( V k + dNi = efjiE 00 - u) H J°°). (120) 

As in the classic Penrose process, the necessary condition for the energy gain is: 

lu h J°° < E™ < 0. (121) 



Thus, in a way fully analogous to the Penrose process for particles (18), one may say that if the 
energy at infinity increases because the black hole absorbed negative-at-infinity energy, then the 
black hole rotation must also slow down by absorbing matter with negative angular momentum. 

Blandford and Znajek [49] made the brilliant discovery that an electromagnetic form of the 
Penrose process may work. In their model, the energy for the jet is extracted from the spin energy 
of the black hole via a torque provided by magnetic field lines that thread the event horizon or 
ergosphere. The estimated luminosity of the jet is given by [178J (although see [303] for higher 
order expressions that apply when a/M ~ 1) 

Lbz = ^ 2 F Blr%{a/Mf , (122) 

where uJp = Vtp{^H — Qf)/Q 2 h is a measure of the effect of the angular velocity of the field Sip 
relative to that of the hole SIr = a/ {r 2 H + a 2 ), B± is the magnetic field normal to the horizon, and 
th is the radius of the event horizon pT| . In this model, the only purpose of the disk is to act as 
the current sheet which continually provides magnetic field to the black hole. This last point led to 
one of the main objections to the Blandford-Znajek model: Ghosh and Abramowicz |109j argued 
on astrophysical grounds that accretion disks simply cannot feed the required fields into the black 
hole. However, recent work by Rothstein and Lovelace |267j has countered this claim and suggested 
that indeed the disk can serve this role. There are also more fundamental reservations with the 
Blandford-Znajek model, some of which are presented in [25211155 ]. Such claims and counter-claims 
were for many years characteristic of the uncertainty in the theory of relativistic jets (see [157] for 
a discussion). However, direct numerical simulations may be helping to clarify the picture, as we 



discuss in Section 11.7 Plus, there is now observational evidence suggesting a possible connection 
between black hole spin and jet power, exactly as predicted by the Blandford-Znajek model [220 , 
although again there are countering claims [93j 



45 



11 Numerical Simulations 



In simulating accretion disks around black holes, there are a number of challenging issues. First, 
there is quite a lot of physics involved: relativistic gravity, hydrodynamics, magnetic fields, and 
radiation being the most fundamental. Then there is the issue that accretion disks are inherently 
multi-dimensional objects. The computational expense of including extra dimensions in a numer- 
ical simulation is not a trivial matter. Simply going from one to two dimensions (still assuming 
axisymmetry for a disk) increases the computational expense by a very large factor (10 2 or more). 
Going to three-dimensions and relaxing all symmetry requirements increases the computational 
expense yet again by a similar factor. Simulations of this size have only become feasible within 
the last decade and still only with a subset of the physics one is interested in and usually with a 
very limited time duration. 

Another hindrance in simulating accretion disks is the very large range of scales that can be 
present. In terms of a grid based code, a disk with a scale height of H/R requiring N z zones to 
resolve in the vertical direction at some radius i?; n , would require something of the order N Z /(H / R) 
zones to cover each factor of R m that is treated in the radial direction. The azimuthal direction in 
a full three-dimensional simulation would require a comparable number of zones to what is used 
in the radial direction. Given that a very large calculation by today's standards is 10 9 - 10 10 zones, 
we can see that treating very thin disks (H/R < 0.01) in a full three-dimensional simulation, 
even with a very modest number of zones in the vertical direction (N z > 50) would take all the 
resources one could muster. However, two-dimensional (axisymmetric) simulations of thin disks 
and three-dimensional simulations of thick, slim, or moderately thin disks have now become quite 
common. 

The approach of numericists in many ways parallels that of theorists, so we structure this 
section much the same as the first part of this Living Review. That is, after a brief introduction 
to numerical techniques, we discuss how the various components of the matter description (a la 
Section[3]) are implemented in numerical simulations. We then review a few special cases illustrating 
how analytic models (Sections [4]-[7]) and numerical simulations can complement one another. We 
finish this numerical section with a few topics of special interest. 

11.1 Numerical techniques 

11.1.1 Computational fluid dynamics codes 

There are many numerical codes available today that include relativistic hydrodynamics or MHD 
that are, or can be, used to simulate accretion disks. A partial list includes: Cosmos++ [18] . 
ECHO [74], HARM [105], and RAISHIN [208] plus the codes developed by Koide et al. [152] . 
De Villiers and Hawley [71], Komissarov [154] . Anton et al. [TS], and Anderson et al. [IB]. This 
represents tremendous growth in the field as prior to 1998 there were only a couple such codes, 
largely derived from the early work of Wilson [318] and later developments by Hawley |1261 1125] . 
Those early codes were only available to a handful of researchers. This was simply a reflection of the 
fact that the computational resources available to most researchers at that time were insufficient 
to make use of such codes, so there was little incentive to develop or acquire them. This has clearly 
changed. 

Today there are primarily two types of numerical schemes being used to simulate relativistic 
accretion disks, differentiated by how they treat discontinuities (shocks) that might arise in the 
flow: the artificial viscosity scheme, still largely based upon formulations developed by Wilson [318J : 
and Godunov-type approaches, using exact or approximate Riemann solvers. Both are based on 
finite difference representations of the equations of general relativistic hydrodynamics, although 
the latter tend to also incorporate finite-volume representations. The artificial viscosity schemes 
have the advantages that they are more straightforward to implement and easier to extend to 
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include additional physics. More importantly, they are computationally less expensive to run than 
Godunov schemes. Furthermore, recent work by Anninos and collaborators [T71 [TS] has shown 
that variants of this scheme can be made as accurate as Godunov schemes even for ultrarelativistic 
flows, which historically was one of the principal weaknesses of the artificial viscosity approach. 
Godunov schemes, on the other hand, appreciate the advantage that they are fully conservative, 
and therefore, potentially more accurate. They also require less tuning since there are no artificial 
viscosity parameters that need to be set for each problem. More thorough reviews of these two 
methods, with clear emphasis on the High-Resolution Shock-Capturing variant of the Godunov 
approach, are provided in the Living Reviews by Marti & Midler [182] and Font ]§T\. Other 
numerical schemes, such as smooth-particle hydrodynamics (SPH) and (pseudo-)spectral methods 
are less well developed for work on relativistic accretion disks. 

11.1.2 Global vs. shearing-box simulations 

Along with settling on a numerical scheme, a decision must also be made whether or not to try to 
treat the disk as a whole or to try to understand it in parts. The latter choice includes "shearing- 
box" simulations, the name coming from the type of boundary conditions one imposes on the 
domain to mimic the shear that would be present in a real disk |121j . The obvious advantage 
of treating the disk in parts is that you circumvent the previously noted problem of the large 
range in scales in the disk by simply ignoring the large scales. Instead one treats a rectangular 
volume generally no larger than a few vertical scale heights on a side and in some cases much 
smaller. In this way, for a moderate number of computational zones, one can get as good, or often 
much better, resolution over the region being simulated than can be achieved in global simulations. 
The obvious disadvantage is that all information about what is happening on larger scales is lost. 
By construction, no structures larger than the box can be captured and the periodic boundaries 
impose some artificial conditions on the simulations, such as no background gradients in pressure 
or density and no flux of material into or out of the box, meaning the surface density remains 
fixed. Perhaps most important in the context of a review on relativity is that the box is treated as 
a local patch within the disk, with a scale smaller than that of the curvature of the metric. Thus 
most general relativistic effects are not, and by construction need not be, included in shearing box 
simulations. Since this is a Living Review in Relativity, we will not dwell much on this class of 
simulations. Still, there are some significant highlights that should be mentioned. 

By far the most extensive work using shearing box simulations has been aimed at better under- 
standing the magneto-rotational instability (MRI). Starting from the earliest "proof-of-concept" 
simulations |121l 1296) , shearing box simulations have been used to demonstrate various properties 
of the saturated state of MRI turbulence [2711 12721 11731 1283j and much about the vertical struc- 
ture of MRI turbulent disks [H] . Shearing box simulations have also proven valuable in studying 
radiation-dominated disks. For example, Turner |307j studied the vertical structure of radiation 
dominated accretion disks and Hirose et al. [129) studied the gas-pressure dominated case, both in- 
cluding radiation using flux-limited diffusion. Hirose and collaborators subsequently used shearing 
box simulations to demonstrate that radiation-dominated disks are thermally stable 11281 , though 
they had long been held not to be [1761 [280 , 246 . As we said, though, there are limits imposed by 
the finite size of the shearing box. For instance, the viscous, Lightman-Eardley instability |176] 
can only be studied through global simulations. 

11.2 Matter description in simulations 

The minimum physics required for a global simulation of an accretion disk are gravity and hydro- 
dynamics (assuming the disk is dense enough for the continuum approximation to hold). Since 
many disks have masses that are small compared to the mass of the central compact object, the 
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self-gravity of the disk can often be ignored. Therefore, in the next three sections, gravity will 
simply mean that of the central black hole. 



11.2.1 Hydrodynamics + gravity 

The first researcher to develop and use numerical algorithms for simulating relativistic accretion 
flows was Wilson 316 , who considered the spherical infall of material with a non-zero specific an- 
gular momentum toward a Kerr black hole using the full metric, although restricted to two spatial 
dimensions. Wilson was able to confirm the additional centrifugal support that the infalling mate- 
rial experienced due to the rotating black hole. For sufficiently high values of angular momentum, 
the material could not immediately accrete onto the black hole, instead forming a fat disk of the 
type described in Section [4] 

In many ways, Wilson's pioneering work was at least a decade ahead of its time and there 
was consequently a lull in activity until Hawley and Smarr collaborated with Wilson to revive his 
work [1261 1125j . As an example of how analytic and numerical work can complement each other, 
it is worth noting that one of the test cases they used for their new two-dimensional relativis- 
tic hydrodynamics code was based on the analytic theory for relativistic thick disks, which had 
been worked out in the time since Wilson's original simulations. Using the analytic theory, they 
constructed a series of disks with different (constant) specific angular momenta, from I < £ ms to 
£ > Cb- The I > £ m b case yielded a static solution as expected and confirmed the ability of their 
code to accurately evolve such a solution in multi-dimensions over a dynamical time. The £ < £ m b 
cases showed greater time variability and illustrated the power of direct numerical simulations to 
extend our understanding of black hole accretion. 



11.2.2 Magnetohydrodynamics + gravity 

Magnetic fields can play many important roles in relativistic accretion disks, from providing lo- 
cal viscous stresses through turbulence that results from the magnetorotational instability (Sec- 
tion 8.2 1, to providing a mechanism for launching and confining jets (Section 10 1. Thus, the 
inclusion of magnetic fields in numerical simulations of relativistic accretion disks is important. 
Although the techniques for doing this were first described in |317j . relatively little work was 
done in this area until fairly recently. Perhaps the two most important contributions so far from 
relativistic MHD simulations of accretion disks have been: 1) elucidating the behavior of MRI 
turbulent disks in the vicinity of black holes, and 2) exploring the many interrelations between 
magnetic fields and relativistic jets. 

In at least one case, the inclusion of MHD + gravity is sufficient to adequately capture the 
dynamics of a real black hole accretion disk. Sgr A*, the black hole at the center of the Milky 
Way galaxy has such an anemically low luminosity |194) that numerical simulations that ignore 
radiation can safely be applied [75] , as has been done by many authors |112l 12271 12111 175] . In most 
other cases, radiative processes must somehow be accounted for. 



11.2.3 Radiation-Magnetohydrodynamics + Gravity 

Probably the most glaring shortcoming of almost all numerical simulations of accretion disks and 
many other phenomena in astrophysics to date is the unrealistic treatment of radiation, which is 
most often simply ignored. This is not due to a lack of appreciation of its importance on the part 
of numericists, but simply a reflection of the fact that there are very few efficient ways to treat 
radiation computationally in multiple dimensions. 

The optically-thin limit is an exception. In this case, the radiative cooling only enters the 
energy and momentum conservation equations as a source term 

V„T^ = fu", (123) 
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where / is the cooling function of the gas (Section 3.4 ) . The subsequent propagation of the radiation 



through the matter can be ignored or post-processed. Even so, including radiative cooling can be 
critically important to capturing the true dynamics of the disk |101l I76j . as shown in Figure 
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In the optically-thick limit, some progress can be made by employing the flux- limited diffusion 
approximation |174j . Global simulations of disks including a flux-limited diffusion treatment of 
radiation (but using pseudo-Newtonian gravity) are now being done by Ohsuga and collabora- 
tors [2321 1231] . A few steps toward the goal of relativistic radiation MHD simulations of black hole 
accretion disks have also been taken in recent years. A method for treating optically thick accre- 
tion using a conservative, Godunov scheme was developed by Farris and collaborators [91] . The 
same basic method has now been used to examine both Bondi [100] and Bondi-Hoyle 326, 266 
accretion. Simulations of accretion disks, though, must await the generalization of this method to 
treat radiation both in the optically thick and thin limits. 

11.2.4 Evolving GRMHD 

In most simulations of accretion disks around black holes, the self-gravity of the disk is ignored. 
In many cases this is justified as the mass of the disk is often much smaller than the mass of the 
black hole. This is also much simpler as it allows one to treat gravity as a background condition, 
either through a Newtonian potential or a relativistic metric (the so-called Cowling approximation 
in relativity). However, there are plausible astrophysical scenarios in which this approximation is 
not valid. Two of the more interesting are: 1) a tidally disrupted neutron star accreting onto a 
stellar-mass black hole; and 2) an overlying stellar envelope accreting onto a nascent black hole 
during the final dying moments of a massive star. Interestingly these scenarios are currently the 
most popular models of gamma-ray bursts 82, 319J, which are the most powerful explosions in the 
Universe since the Big Bang and are observed at a rate of about one per day. 

Treating a self-gravitating disk in general relativity requires a code that can simultaneously 
evolve the spacetime metric (by solving the Einstein field equations) and the matter, magnetic, 
and radiation fields. Codes capable of doing this have very recently become available, divided into 
the following "flavors": self-gravity + hydro self-gravity + MHD [STJ HH EH] ; and self- 
gravity + radiation MHD (91]. More detailed discussion of these methods and their application 
to problems in astrophysics is provided in the Living Review by Font [97j . 

The importance of evolving the gravitational field in the case of black hole accretion is well 
illustrated in the case of the runaway instability, which we described in Section |8.1.2| That in- 
stability has now been simulated directly, including evolving the spacetime metric [2091 12101 1160) . 
These studies have confirmed earlier results that the runaway instability does not develop. 

Another application of these types of codes that is relevant to our review is the reexamination of 



the Papaloizou-Pringle instability (discussed in Section 8.1.1 1. Through this instability, a massive, 
self-gravitating Polish doughnut orbiting around a black hole could become a strong emitter of 
large amplitude, quasiperiodic gravitational waves j!49j . 

11.3 Polish doughnuts (thick) disks in simulations 

As we already mentioned, it is computationally less expensive to simulate thick disks than thin, so 
it comes as no surprise that most of the numerical work for nearly four decades, starting with the 
work of Wilson |316j . has been on thick disks. In this section we touch on a few highlights from 
this area of research 

Some of this work has focused on confirming predictions from analytic theory. Igumenshchev 
and Beloborodov [137] used two-dimensional relativistic hydrodynamic simulations to confirm that: 
(i) The structure of the innermost disk region strongly depends on the black hole spin, and (ii) 
the mass accretion rate is proportional to the size of the energy gap between the inner edge of the 



disk and the cusp (see Figure 15), as predicted by Kozlowsk and collaborators 
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Figure 14: Pseudo-color plots of log(T) with contours of log p from 4 different general relativistic 
MHD simulations. The simulations all begin with the same initial conditions, but have different 
energy conservation and cooling treatments: The upper-left panel conserves internal energy and 
ignores cooling; the upper-right panel conserves internal energy and includes cooling; the lower-left 
panel conserves total energy and ignores cooling; the lower-right panel conserves total energy and 
includes cooling. The very different end states illustrate the importance of properly capturing 
thermodynamic processes. Image reproduced by permission from |101j . copyright by AAS. 
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Figure 15: Time-average mass accretion rate rh = M^/L^d as a function of the energy gap A<I> 
for models with a = (circles), a/M = 0.9 (squares), and a/M = —0.9 (triangles). The bars show 
the variability of m. The lines represent the predicted dependencies rn oc 

(A$)W(7 9 -i) where 

7 9 = 4/3 is the adiabatic index. Image reproduced by permission from |137j . copyright by RAS. 
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Hawley (1 161 II IT] studied the nonlinear evolution of the non-axisymmetric Papaloizou-Pringle 



instability (Section 8.1.1). He discovered that for radially slender tori, the principal mode of 



the PPI saturates with the formation of ellipsoidal overdensity regions he termed "planets" (see 



Figure 16 for an example). In cases where several wavenumbers were unstable, multiple planets 
would begin to form but then nonlinear mode-mode coupling would cause them to merge. For 
wide tori, particularly when the tori extended beyond their Roche limits, overflowed their cusps, 
and accreted into the hole, the PPI saturated at low amplitude or was prevented from growing 
altogether. The initial work only considered Schwarzschild black holes, but all of the results were 
later qualitatively confirmed for Kerr black holes |70) . 

Another interesting note is that most of the global numerical simulations of MPJ turbulent 
disks have used a Polish doughnut as the starting condition. This was first done by De Villiers and 
Hawley [55J (animations can be viewed at [6T)Illl5| ) and later by a number of other groups [1051118) . 
As expected from linear analysis and earlier Newtonian simulations, the Maxwell stresses produced 
by the MRI-driven turbulence naturally force the final distribution of the angular momentum to 



be nearly Keplerian (see Figure 17), independent of the initial matter distribution. 



Although the magnetic field plays a crucial role through the action of the MRI, its amplitude 



saturates at a relatively weak level, always remaining subthermal in the disk (see Figure 18). On 
the other hand, simulations show that magnetic fields do tend to dominate regions outside the 
disk, including in the hot, magnetized corona that sandwiches the disk and in the evacuated, 
highly magnetized funnel region. The funnel region and the boundary between the funnel and the 



corona are where jets and outflows arc gencrically observed in these simulations (Section 11.7). 

Despite the increased complexity of the disk when MRI-driven turbulence is at play, many 
features and properties remain strikingly similar to what is revealed in hydrodynamic-only stud- 
ies. For instance, a commonly seen feature in global, non-radiative MHD simulations of black 
hole accretion disks is an "inner torus" [120] . The inner torus is predominantly a hydrodynamic 
structure, as it is largely gas pressure supported. Generically, this inner torus consists of a se- 
quence of equidensity and equipressure surfaces with a pressure maximum at ~ 10 r^, a cusp at 
~ 3t"g, and a roughly parabolic, evacuated, though magnetized, funnel along the rotation axis. 
This structure is remarkably similar to the Polish doughnut model, as shown in Figure [19] Thus, 
many hydrodynamic torus results may retain relevance even in MRI turbulent disks [253] . 

11.4 Novikov-Thorne (thin) disks in simulations 

Early numerical simulations of black hole accretion disks nearly all focused on thick disks. This was 
mostly for computational convenience as thicker disks require fewer resources than thin ones. In 
recent years, though, the resources have become available to start testing thinner disks. This has 



enabled researchers to begin rigorously testing the Novikov-Thorne model (Section 5.3), especially 
the assumption that the internal torque of the disk vanishes at the ISCO. This is important as 
some researchers have argued that magnetic fields might nullify this hypothesis by maintaining 
stresses inside the ISCO [T64 ] IT04 ] [25] . 

Initial global MHD simulations did show some variations from the Novikov-Thorne model at 
the ISCO at the level of - 10% [T23~l [2~2~6] . while others [M5] showed much smaller deviations < 3%. 
Penna and collaborators |240j discussed possible reasons for the discrepancies, and suggested that 
the results are, in fact, in better agreement than they may have initially appeared. Their conclusion 



is that, in fact, the Novikov-Thorne model matches numerical results quite well (see Figure 20 plus 
animations can be viewed at [239] ) . Any remaining differences are small enough to have negligible 
effect on common applications of the Novikov-Thorne model, such as measuring black hole spin 



via spectral fitting [TBB][327j (Section 12.1) 



We mention that since present day relativistic numerical simulations do not treat the radiation 
in geometrically thin, optically thick disks, the simulations of Novikov-Thorne disks have so far 
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Figure 16: Equatorial slice through hydrodynamic tori at saturation of the Papaloizou Pringlc 
instability showing formation of significant non-axisymmetric (m = 1) overdensity clumps. The 
density contours are linearly spaced between /0§gix and 0.0. This figure represents models A3p 
(left) and B3r (right) of (70]. 




Figure 17: Specific angular momentum £ as a function of radius at t = (thin line) and at 
t = 10.0 orbits (thick line). The individual plots are labeled by model. In each case the Keplerian 
distribution for a test particle, ^Kep, is shown as a dashed line. Image reproduced by permission 
from [72!, copyright by A AS. 
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Figure 18: Color contours of the ratio of azimuthally averaged magnetic to gas pressure, -P m ag/-Pgas- 
The scale is logarithmic and is the same for all panels; the color maps saturate in the axial funnel. 
The body of the accretion disk is identified with overlaid density contours at 1CP 2 , 10~ 15 , 10 _1 , 
and 10 -0 5 of p max (i = 0). The individual plots are labeled by model. In all cases, the magnetic 
pressure is low (P m ag/-Pgas <C 1) in the disk, comparable to gas pressure (-P mag /-Pgas ~ 1) in the 
corona above and below the disk, and high (P m ag/-Pgas 3> 1) in the funnel region. Image reproduced 
by permission from [72], copyright by A AS. 
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Figure 19: On the left, equidensity contours calculated from an analytic Polish doughnut. On 
the right, equidensity contours from a numerical MHD simulation (model 90h from [H]). Note, 
though, that the contours on the left are linearly spaced, while those on the right are logarithmically 
spaced. Thus, the gradients represented on the left are shallower than those on the right. Image 
reproduced by permission from [253] . copyright by ESO. 
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Figure 20: Left: Time-averaged rest mass density in the r-z plane for four GRMHD simulations 
with a = and various disk thicknesses. The dashed vertical line marks the ISCO. The disk 
opening angle, h = H/r, and effective Shakura-Sunyaev viscosity, a, are reported in each panel. 
The top three panels have h <C a and the inner edge of the disk is located outside the ISCO. 
The bottom panel has /i> a and the density increases monotonically down to the event horizon. 
Figure from [241] . Right: Various fluxes as functions of radius for a numerical Novikov-Thorne 
disk simulation. Top: Mass accretion rate. Second panel: Accreted specific angular momentum. 
Solid line is simulation data; dashed line gives Novikov-Thorne solution; dotted line is ISCO value. 
Note that the specific angular momentum does not drop significantly inside the ISCO, unlike for 



thick disks, such as in Figure 17 Third panel: The "nominal" efficiency, which is the total loss of 
specific energy from the fluid. Bottom panel: Specific magnetic flux. The near constancy of this 
quantity inside the ISCO is an indication that magnetic stresses are not significant in this region. 
Image reproduced by permission from [240J . 
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employed an ad hoc cooling prescription [2781 [226! . This prescription conveniently makes the 
same assumption that the Novikov-Thorne model does: that all energy dissipated as heat in the 
disk is radiated away locally on roughly the orbital timescale. This is probably reasonable for an 
appropriate range of mass accretion rates, though it will be good to test this assumption with 
future global radiation-MHD simulations. 



11.5 ADAFs in simulations 

A lot of recent simulation work has been focused on exploring more ADAF-like flows under the 
action of realistic MM turbulence [222 , 322 , 32T)|. To some extent, this is an extension of the Polish- 
doughnut simulations of the last decade, yet goes beyond it in at least two important respects: 1) 
the simulations cover a significantly larger spatial range (a few hundreds of GM/c 2 versus a few 
tens); and 2) the simulations explore much longer temporal evolution (hundreds of thousands of 
GM/c 3 versus tens of thousands). The result of this is that the simulations are able to explore 
steady-state accretion out to much larger radii (~ 100 GM/c 2 as opposed to ~ 10 GM/c 2 ). More 
can be expected from this work in the coming years. 



11.6 Oscillations in simulations 

Hydrodynamic simulations have also been used extensively to study the natural oscillation modes 



of relativistic disks, particularly as they might relate to QPOs (Section 12.4). Rezzolla and col- 
laborators |263[ I264j identified a mode they referred to as a p-mode that occurred in a near 3:2 
ratio with the radial epicyclic mode, which was subsequently confirmed through numerical simu- 
lations |3251 1324] (animations can be viewed at |262J). The 3:2 ratio of these modes is important 
as the highest-frequency QPOs in black hole low-mass X-ray binaries are observed to occur in this 
ratio. Later, Blaes and collaborators [IS] identified a different pair of modes, the vertical epicyclic 
and axisymmetric breathing mode, that have a near 3:2 ratio for a broader range of parameters. 
Similar to |325l 1324) , they also demonstrated that these modes could be identified in numerical 
simulations. 

A number of MHD simulations have also been performed which claim to observe QPOs. 
Kato |145) performed a global MRI simulation in a pseudo-Newtonian potential and claimed to see 
QPOs in the power spectrum in a measure of the luminosity associated with radial infall. Machida 
and Matsumoto [179] have reported the formation of a one-armed spiral within the inner torus of 
a global MRI simulation, and suggested that this could be responsible for the low frequency QPO 
in the hard state. Schnittman, Krolik, and Hawley [274J have found tentative evidence for high 
frequency QPO features in light curves generated by coupling a general relativistic ray tracing code 



to a global GRMHD simulation with an assumed emission model. Henisey and collaborators [127] 
found similar tentative evidence in a sample of tilted disk simulations, possibly confirming earlier 
suggestions that disk tilt may be an important mechanism for driving high frequency QPOs in 
black hole accretion disks |1431 |9"5] . Most recently Dolence and collaborators [75] report evidence 
for near- infrared and X-ray QPOs in numerical simulations of Sgr A*. 

Aside from these few interesting examples, though, no robust QPO signal has generically 
emerged in MHD simulations 22, 261] . It is unclear what the implications of this are. It may be 
that some missing physics, such as radiation transport, plays a fundamental role in exciting QPOs. 
This is an important open problem in numerical simulations of black hole accretion disks. 



11.7 Jets in simulations 

In recent years several numerical simulations have demonstrated the generation of jets self-consistently 
from simulations of disks. A few researchers have claimed to produce jetted outflows from purely 
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hydrodynamic interactions. Nobuta and Hanawa |228j , for instance, were able to achieve jetted 
outflows driven by shock waves when infalling gas with high specific angular momentum collided 
with the centrifugal barrier. However, such simulations usually require very special starting condi- 
tions and the jets tend to be transient features of the flow; therefore, it is unlikely these are related 
to the well-collimated expansive jets observed in many black hole systems j^j 

Among MHD simulations, some distinction should be drawn between those that impose large 
scale magnetic fields that extend beyond the domain of the simulation and those that, at least 
initially, impose a self-contained magnetic field. In the first case, the disk is often treated merely as 
a boundary condition for the evolution of large-scale magnetic field. Shibata and Uchida 284, 285 
used this technique to show that jets could be driven both by a pinching of the fields due to radial 
inflow through the disk and by the j x B force caused when twisted fields unwind. 

Since 1999, many MHD disk simulations, usually starting with locally confined, weak magnetic 
fields, have shown some sort of magnetically driven bipolar outflow [1531 11461 1731 11241 11911 [189 . 
Since these jets are propagating at the Alfven speed in a magnetically dominated region (/3 m <C 
1), they proceed supersonically relative to the gas. However, some authors claim the outflow is 
only a transient state |146| and most do not achieve the very high Lorentz factors observed in 
astrophysical jets [73] • It is also unclear whether or not the jets in these simulations are connected 



with the Blandford-Znajek process [73] (Section 10 1, although see |188j . McKinney [189] and 
Komissarov |155j suggest that part of the problem is that the jets arising from MHD simulations 
of disks are often contaminated by the application of numerical floors on density and energy and 
through numerical diffusion, both of which can lead to anomalous mass loading of the jet. In their 
own carefully constructed simulations, they independently found that the magnetic structures 
exhibited by their simulated jets are, in fact, consistent with the expectations of the Blandford- 



Znajek model, as illustrated in Figure 21 (animations can be viewed at [187] V 

Tchekhovskoy and collaborators [303] have used GRMHD simulations of black hole accretion 
disks plus jets to investigate possible explanations for the observed radio loud/quiet dichotomy 
m AGN0 For a black hole surrounded by a thin disk, the Blandford-Znajek mechanism predicts 
the luminosity of the jet should scale as Lbz °t Q 2 h- This limits the range of power expected for 
realistic AGN spins to a factor of a few tens at most - too small to explain the observed differences. 
Thicker disks, on the other hand, produce jets whose luminosity can scale as Lbz °c Sl A H or even 
Q% [303 , providing sufficient range to perhaps explain the dichotomy. 

There has also been some work recently trying to understand a different jet dichotomy - one 
that is observed in black hole low-mass X-ray binaries (LMXBs). These exhibit radio emission 
(associated with jets) whose properties change with the observed X-ray spectral and, to a less well 
determined extent, temporal properties of the accretion disk [551 [SI]. Briefly, compact, steady jets 
are observed in the Low/Hard state, whereas jets appear absent in the High/Soft state. Fragile 
and collaborators [102] used numerical simulations to rule out disk scale height as the controlling 
factor, suggesting instead that perhaps the jets are intimately connected with the corona or failed 
MHD wind. Alternatively, it could be that magnetic field topology is the key factor 1361 1190] , 

11.8 Highly magnetized accretion in simulations 

Recently, work has begun to focus on highly magnetized disk configurations, for which some modes 
of the MRI may be suppressed. One motivation is that these might provide an alternate explanation 



for the Low/Hard state [35 (Section 12.3) 



13 For the benefit of readers who are unfamiliar with the phenomenology of relativistic jets we mention that many 
jets demonstrate incredibly consistent collimation, even along hundreds or thousands of kiloparsecs (e.g. |114| ). 
From this we can infer that such jets have maintained their orientation and outflow for millions of years. 

14 Radio loud AGN are ~ 10 3 - 10 4 times brighter in radio than radio quiet AGN of comparable optical luminosi- 
ties [2891 . 
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Figure 21: On the left is a schematic diagram of the Blandford-Znajek mechanism [15] for an 
assumed parabolic field distribution. On the right is the result of a numerical simulation from [158 
showing a very similar structure. Images reproduced by permission; copyright RAS. 
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One way a highly magnetized state could come about is as a result of a thermal instability in an 
initially hot, thick, weakly-magnetized disk [1791 1101] . Machida and collaborators |179j simulated 
this process for an optically thin disk, assuming bremsstrahlung cooling and pseudo-Newtonian 
gravity. They found that, indeed, the densest inner regions of the disk collapse down to a cool, 
thin, magnetically-supported structure. Fragile and Meier |101) extended these results by including 
more cooling processes and using a general relativistic MHD code. 

Another way to achieve a highly magnetized state is to have the disk "drag" the magnetic field in 
from some distant region [121 1291). If the field has a consistent net flux, then doing so will necessarily 
increase the strength of the field near the hole due simply to the smaller area through which the 
flux must thread. Recent numerical simulations have demonstrated this effect, first in the pseudo- 



Newtonian limit [135] (Figure 22 1, and more recently in relativistic simulations |304[ 1192] , These 
examples are interesting because they have led directly to a phenomenon known as a "magnetically 
arrested" accretion state |218j . where the accumulation of magnetic field near the black hole is 
sufficient to temporarily halt the inflow of matter. This state has been shown to produce jet 
efficiencies rj = (M - E)/(M) > 1 [304], where 



M(t) = - J y^)pu r d0d4> (124) 

is the mass accretion rate and 

E{t) = { y/=gT[d9d<t) (125) 



is the energy flux, both taken at the black hole event horizon rn, and the angle brackets indicate 
a time-averaged quantity. This can only happen if at least some of the energy powering the jet is 
being extracted from the rotational energy of the black hole itself! 
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Figure 22: Top: Distributions of density in the meridional plane at different simulation times, 
showing a magnetically arrested state (left) and a non-arrested state (right). Bottom: Snapshot 
of magnetic field lines at the same simulation times. Image reproduced by permission from |135| . 
copyright by A AS. 
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12 Selected Astrophysical Applications 



12.1 Measurements of black-hole mass and spin 

Astrophysical black holes are not charged, and thus are characterized only by their mass and spin. 
Measurements of black hole mass are generally straightforward, requiring only the observation of 
an orbital companion and an application of Kepler's laws. Current mass estimates for stellar-mass 
black holes (in particular microquasars) are reviewed by McClintock & Remillard [184J and for 
supermassive black holes by Kormendy & Richstone |159j (see also j!96j for somewhat more recent 
data). 

Nevertheless, there remain several fundamental questions connected with the mass of black 
holes, and some of them are directly connected to accretion disk theory. One of them is the 
question of ultra- luminous X-ray sources (ULXs) [611 H81J . ULXs are powerful X-ray sources 
located outside of galactic nuclei, which have luminosities in excess of the Eddington limit for 
M = 10 2 Mq, assuming isotropic emission. The huge luminosities of ULXs lead some to conclude 
that they are accreting intermediate-mass black holes, with masses Mulx > 10 2 Mq (e.g. [1651 
198 ). Others think that they are stellar mass black holes with Mulx ~ 10 M Q , either exhibiting 



beamed emission 108 or surrounded by disks that are somehow able to produce highly super- 
Eddington luminosities (e.g. [147j|36]). At this time, at least one ULX (ESO 243-49/HLX-l) has 
been convincingly demonstrated to have an intermediate mass (~ 500 Mq) [501 157] . 

As difficult as it is to nail down the mass on some of these systems, it is even more difficult to 
measure black hole spin, even though it plays a direct, and important, role in accretion disk physics. 
One obvious example of the role of spin is the dependence of r ms , the coordinate radius of the ISCO, 
on spin. In the symmetry plane of the black hole r ms — Qtq for a* = a/M = (non-rotating), 
Itq for a* = 1 (maximal prograde rotation) , and 9 tq for a* = — 1 (maximal retrograde rotation) 



(Section 2.5). It is believed that the inner edge of the accretion disk will be similarly affected. In 
addition, for rapidly rotating black holes, the Blandford-Znajek mechanism, and similar processes 
which depend on spin, may account for a fair share of the global energetics, comparable to that 
of accretion itself (see e.g. |155| I189j and references therein). Therefore, measuring the spin of 
accreting black holes is integrally tied up in understanding black hole accretion generally. 

Four methods for determining black hole spin have been proposed in the literature. With 
references to some of the earliest results, they are: 1) fitting the continuum spectra of microquasars 
observed in the thermally dominant state using disk emission models [5S1 11851 11971 1277] ; 2) fitting 
observed relativistically broadened iron line profiles with theoretical models |142l 13151 12011 12021 
1260] : 3) matching observed QPO frequencies to those predicted by theoretical models [62" 1 [TT | |259[ 
306 ; and 4) analyzing the "shadow" a black hole makes on the surface of an accretion disk |301) . 
The first three methods have been the most commonly applied to date. Although there have been 
some glaring discrepancies in the spin estimates published to date (e.g. one group claiming that 
Cyg X-l has a near-zero spin, a* = 0.05 ±0.01 [200] . and then later claiming it has a near-maximal 
spin, a* > 0.95 [89), there appears to be a settling of values in recent years and a growing 
confidence in the methods, particularly the continuum fitting. There are still some concerns, 
however. For the continuum fitting method, the main issue is that the inclination of the X-ray 
emitting region must be measured by some independent means. This is because the effect of the 
inclination on the spectrum is degenerate with the effect of spin [175] , so both can not be accounted 
for within the continuum fitting method. In cases where such an independent measure is available 
(e.g. |292j ). the continuum fitting method appears robust. For the relativistically-broadened iron 
line method, there are difficulties in properly estimating the extent of the "red" wing, which is 
most directly related to the spin of the black hole, and in modeling the hard X-ray source photons 
and the disk ionization, both of which strongly affect the reflection spectrum. The reviews by 
Remillard, McClintock, and collaborators [2581 1186j give more complete introductions to the topic 
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of measuring black hole spin, with emphasis mainly on the continuum fitting method. 



12.2 Black hole vs. neutron star accretion disks 

Little of what we have said so far has depended on whether the central compact object is a black 
hole or neutron star, provided only that the neutron star is compact enough to lie inside the inner 
radius of the disk r- m . In this case, its presence will not be noticed by the disk except through 
its gravity, which will be practically the same as for a black hole (an exception would be if the 
neutron star is strongly magnetized |113j ). However, this does not mean that accreting black hole 
and neutron star sources will be indistinguishable, as we have not yet fully addressed the question 
of what happens to energy advected past r- ln . For optically thick, geometrically thin Shakura- 



Sunyaev disks (Section 5.3), a significant fraction of the gravitational energy liberated by advection 
is radiated by the gas prior to it passing through r m . Thus, the total luminosity of thin disks will 
not depend sensitively on the nature of the central object. However, this is not the case for the 
ADAF solution (Section [7]), for which much of the thermal energy gained by the gas from accretion 
is carried all the way in to the central object. Narayan and his collaborators [2251 T215, 195, 1107) 
have convincingly argued that this may allow observers to distinguish between black hole and 
neutron star sources. 

The key is that, for black hole sources, advection through the event horizon allows the excess 
thermal energy to be effectively absorbed without ever radiating. For neutron star sources, on the 
other hand, the presence of a hard surface ensures that the excess energy of accretion is released 
upon impact and must be radiated to infinity. This implies that for systems in the ADAF state, a 
black hole source should be significantly less luminous that a neutron star one with the same mass 
accretion rate |225j . Perhaps a more important point is that the range of luminosities should be 
wider for a black hole source than for a neutron star one |215j . This is because, while the luminosity 
goes as L oc to for all neutron star states and for black holes in a high accretion state, it goes as 
L cx to 2 for black holes in the ADAF state for which m < 10~ 2 — 10 . Furthermore, since the 
luminosity is also proportional to the mass of the central object L oc M, at the highest accretion 
rates a black hole source should be more luminous than a neutron star one due to its higher mass 



(Figure 23). Another key point to this argument is that neutron stars can independently and 
reliably be confirmed if they display type I bursts, which are thermonuclear flashes occurring in 
material accumulated on the surface of the neutron star [141] . Thus one can compare known 
neutron star sources against suspected black hole sources. This has now been done in a number of 



otherwise similar sources and Narayan's expectations have indeed been confirmed 215, 195, 107 



(a recent example is shown in Figure 23 ) . This provides compelling observational evidence for the 
existence of black hole event horizons, although this falls short of being a proof [T^]. This topic is 
discussed further in the review article by Narayan and McClintock |219j . 



12.3 Black-hole accretion disk spectral states 

Black hole accretion disks, particularly in X-ray binaries, exhibit complex spectra composed of both 
thermal and nonthermal components. During outbursts, the relative strengths of these components 
change frequently in concert with changes in luminosity and the characteristics of the radio features 
(i.e., jets). Astronomers have developed a set of empirical spectral classification states to broadly 
characterize these observations. These spectral states likely reveal important information about 
the underlying physical state of the system; therefore it is worth summarizing the states and their 
observed properties here. 

Probably the easiest state to connect with a theoretical model is the "High/Soft" (HS) or 
"Thermally Dominant" state. As the name implies, the spectrum in this state is dominated by 
the thermal component. This state is best explained as ~ 1 keV thermal emission from a multi- 
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Figure 23: Left: Luminosity as a function of accretion rate for neutron star and black hole sources, 
illustrating that a wider range of luminosities are expected for black holes. Image reproduced by 
permission from [215] , copyright by AAS. Right: Recent data showing that neutron star sources 
(open symbols) are systematically more luminous than black hole sources (filled symbols) in anal- 
ogous spectral states. Image reproduced by permission from |171) . copyright by Elsevier. 



65 



temperature accretion disk, as predicted by the Shakura-Sunyaev (thin) disk model (Section 5.3 1. 
However, most sources spend the majority of their lifetimes in the "Low/Hard" (LH) or even "qui- 
escent" state. The quiescent state is characterized by exceptionally low luminosity and a hard, 
non-thermal spectrum (photon index T = 1.5 — 2.1). As the luminosity increases the sources 
usually enters the Hard state. Here the 2-10 keV intensity is still comparatively low and the 
spectrum is still nonthermal. This spectrum is best fit with a powerlaw of photon index T ~ 1.7 
(220 keV). In this state, the thermal component is either not detected or appears much cooler, 
indicating the thin disk may truncate further out than in the Thermally Dominant state, although 
see |199| for claims that the thin disk extends all the way to the ISCO even in the Hard state. 
Observations suggest the region interior to the thin disk may be filled with a hot (presumably 
thick), optically-thin plasma, which accounts for the nonthermal part of the spectrum. This is 
the picture suggested by the "truncated disk model" [55] HEl HI] , which pictures the Hard state 
as a truncated thin (Shakura-Sunyaev) disk adjoined with an inner thick (ADAF-like) flow. This 
model has shown tremendous phenomenological success [75], although other models for the Hard 
state still abound [323], 11011 1136] . The Hard state is also linked with observations of a persistent 
radio jet that is not seen in other states (see Section 10 1. The final spectral state, which is referred 
to as the "Very High" (VH) or "Steep Power Law" state, is characterized by the appearance of 
high-frequency QPOs (Section 12.4) and the presence of both disk and powerlaw components, each 
of which contributes substantial luminosity. In this state, the powerlaw component is observed to 
be steep (r ~ 2.5), giving the state its name. This state is sometimes associated with intermittent 
jets. 

These states, their distinguishing observational properties, and a sampling of observations from 
various black hole X-ray binaries is presented in the review by McClintock and Remillard 184 . 



12.4 Quasi-Periodic Oscillations (QPOs) 

Einstein's general theory of gravity has never been tested in its strong field limit, characteristic 
of the region very near black holes (or neutron stars), i.e., a few gravitational radii away from 
these sources. Soon VLBI measurements may be able to resolve these scales for the supermassive 
black hole at the center of our galaxy. However, for most sources, resolution in time seems to be 
a more practical approach. To zeroth order, the light curves from accreting black holes vary in 
a chaotic manner, resembling a loud noise. However, Fourier analysis of the light curve reveals 
stinkingly regular patterns buried in the noise. For Galactic black hole and neutron star sources, 
these quasi-periodic oscillations (QPOs) have frequencies of a few hundred Hz. 

Frequencies in the range 100-1000 Hz formally correspond to orbital frequencies a few gravita- 
tional radii away from a stellar-mass object. The focus on orbital frequencies is further motivated 
by the stability of observed QPO frequencies over very long periods of time. For example, QPOs 
of 300 Hz and 450 Hz were observed from the microquasar GRO J1655-40 during its 1996 out- 
burst and again almost nine years later during its 2005 outburst. This strongly suggests that the 
frequencies cannot depend on quantities such as magnetic field, density, temperature, or accretion 
rate, as these all vary greatly in time. The only parameters of a black hole accretion system that 
do not vary over a nine year period are the mass and the spin of the central black hole. Thus, the 
oscillation frequencies must only depend on these two parameters, and only frequencies connected 
to orbital motion have the property that they depend only on mass and spin. Thus, the possible 
frequencies are: the Keplerian frequency, the two epicyclic frequencies (as originally suggested by 
Kluzniak and Abramowicz |llj). the Lense-Thirring frequency (as originally suggested by Stella 
and Vietri [294] ). and their combinations (e.g., |145U143j ). 

In several microquasars the detected high-frequency QPOs come not as single oscillations, but as 
part of a pair. Furthermore, Abramowicz and Kluzniak |llj noticed that they are commensurable, 
being most often in a 2/3 ratio, as shown in Table [2] This suggests that a resonance may be at 
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Table 2: Frequency ratio of the "twin peak" QPOs in all four microquasars where they have been 
detected. 



Microquasar Frequency ratio 



GRO J1655-40 


300/450 = 


0.66 


XTE 1550-564 


184/276 = 


0.66 


H 1743-322 


166/240 = 


0.69 


GRS 1915+105 


113/168 = 


0.67 



work. Twin peak QPOs in the kilohertz range have been also detected from binaries containing 
accreting neutron stars. These neutron star QPOs show a similar, though less obvious, 2/3 ratio. 

It was realized j!50j that the behavior of the observed QPO frequencies and amplitudes in 
both neutron star binaries and microquasars is typical for a certain type of non-linear resonance. 
Indeed, it may be observed when a properly tuned spring is attached to a pendulum with a 
properly chosen length. Such a system oscillates in two "modes": the pendulum mode and the 
spring modej^] Mostly due to efforts of Rebusco and Horak [2561 1130L 11311 1132] . a mathematical 
resonance model was developed to describe an arbitrary system oscillating in two modes near a 
2/3 non-linear resonance. The model's predictions for the frequency and amplitude behaviors are 
strikingly similar to the ones observed in X-ray binary QPOs, suggesting that they may indeed 
be explained as a non-linear resonance of two modes of oscillation. The model does not, however, 
explain what these modes are, how they are excited, nor what energy reservoir they tap. Only 
when these questions are answered satisfactorily could one say that the QPO puzzle is solved P'] 

As a final note, a crucial discovery by Barret and collaborators [32j [33l [34] concerning the 
behavior of the quality factor of twin peak QPOs proves that they are disk oscillations and cannot 
be explained by kinematic (Doppler) effects due to the presence "hot spots" on the accretion disk 
surface. These effects are, however, important in modulating the QPO's signal [55]. It is not clear, 



though, whether the oscillations are explained by the discoseismic modes discussed in Section 9.2 
(see, e.g., |313j for references). 



12.5 The case of Sgr A* 

As already mentioned in Section [2. there is a very good chance that the first direct evidence for 
a black hole event horizon will come from Sgr A*, the compact, supermassive object at the center 
of the Milky Way. There have already been a number of strong, indirect arguments in favor of the 
black hole nature of Sgr A* [55] [53], but no direct evidence yet. Sgr A* is also of interest because 
it represents a unique case of black hole accretion, having by far the lowest (scaled) mass accretion 
rate and radiative efficiency of any known source. We anticipate greatly expanding this section in 
the near future as new results become available. 



13 Concluding Remarks 

Since the early 1970s, the study of black hole accretion disks has yielded to remarkable successes. 
And yet, as in most fields of research, each step forward has been met by new questions. In this 

15 It is quite remarkable that the mathematical theory describing resonances makes some very general predictions 
about the behavior of frequencies and amplitudes, even if the specific physical properties of the oscillating objects 
are not known. To some extent it is possible to accurately describe how things oscillate without even knowing what 
is oscillating. 

16 The mathematical foundations of the resonance model are described and discussed in a special issue of the 
Astronomische Nachrichten [2]. 
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article, we have tried to give a tour of some of the successes, such as the many disk models (thick, 
thin, slim, ADAF, . . . ) that have given a firm foundation on which to work, the many studies 
of disk instabilities and oscillations that help us to understand the ways in which real disks can 
deviate from the simplistic models, and the numerical simulations that come as close as possible to 
an experimental test-bed for black hole accretion. We have also tried to indicate what we believe 
are some of the most pressing challenges of the day, including matching our theoretical knowledge 
to actual observed phenomena such as black hole spectral states, quasi-periodic oscillations, and 
relativistic jets. There are also observational challenges to find direct evidence of black hole event 
horizons and definitively constrain black hole spins. We hope, as we continue to update this Living 
Review, to be able to report on future discoveries in these areas, just as we expect to report new 
puzzles we have yet to encounter. 

Clearly our tour has been incomplete. For instance, despite their prominent role in nature, e.g., 
as an AGN feedback mechanism, outflows are not accounted for in any of the four main accretion 
models we presented, as the models all assume that the accretion rate is constant with radius. In 
reality, outflows may be triggered by any of three mechanisms: thermal, radiative, or centrifugal. 
Thermal winds are expected to result from heating of the outer regions of an accretion disk by its 
hot inner region. Radiative winds are driven by radiative flux acting on line opacities. Centrifugal 
acceleration of particles can take place along magnetic field lines which are sufficiently inclined to 
the disk plane. 

The ADIOS (adiabatic inflow-outflow solution) [37] is a generalization of the ADAF solution 
that actually has much of the mass, energy, and angular momentum of the accretion "disk" carried 
away in the form of winds, rather than being advected into the black hole as in a normal ADAF. 
However, the argument behind the ADIOS model is flawed [Hj. Blandford and Begelman [37] 
claim that black hole accretion flows with small radiative efficiency must necessarily experience 
strong outflows because the matter all has a positive Bernoulli constant. Yet a positive Bernoulli 
constant is only a necessary, but not sufficient, condition for outflows. For example, the classical 
Bondi accretion solution has B > everywhere and yet experiences no outflows. Furthermore, low 
efficiency accretion flows do not really have a positive Bernoulli constant everywhere, as boundary 
conditions (ignored in [47 ) will impose some regions of negative Bernoulli constant. A more recent 
look at the ADIOS solution is presented in |37| . 

There are several other analytic and semi-analytic models of accretion disks. Some are closely 
related to the models we have already discussed. For example, the CDAF (convection-dominated 
accretion flow) [2171 1221) is another variant on the ADAF, in which long-wavelength convective 
instabilities transport angular momentum inward and energy outward. Other models relax some 
of the standard assumptions about accretion disks. For example, Bardeen and Petterson [3D] and 
others [1671 12731 1177j relaxed the assumption that disks are axisymmetric by considering tilted 
accretion disks, acted on by the Lense-Thirring precession of the central (rotating) black hole. 
Then there is the exact solution for stationary, axisymmetric non-circular, accretion flows found 
by Kluzniak and Kita |151j . 

On the numerical side, too, there are many interesting accretion configurations that have been 
identified, but are not included in this review. Some examples include quasi-spherical (low angular 
momentum) accretion flows 250, 249] and convection-dominated disks [138] , 

For those who wish for more details or a different perspective, we can recommend several 
excellent text books and review articles devoted, partially or fully, to black hole accretion disks. 
We recall here some of the most often quoted. The oldest, but still very useful and informative, 
is the classic review by Pringle [248 . The most authoritative text book on accretion is Accretion 
Power in Astrophysics by Frank, King and Raine [103] . Two monographs devoted to black hole 
accretion disks are: Black-Hole Accretion Disks by Kato, Fukue and Mineshige [144] . and Theory 
of Black Hole Accretion Disks by Abramowicz, Bjornsson, and Pringle [3J. Lasota |170j also wrote 
an excellent non-technical Scientific American article on black hole accretion in microquasars. 
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Finally, there is a nice series of lecture notes by Ogilvie available on the web |230j . 
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